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Abstract 



Lyapunov exponents, important invariants of a complicated dynamical pro- 
cess, can be difficult to determine from experimental data. In particular, when 
using embedding theory to build chaotic attractors in a reconstruction space, 
extra "spurious" Lyapunov exponents can arise that are not Lyapunov expo- 
nents of the original system. By studying the local linearization matrices that 
are key to a popular method for computing Lyapunov exponents, we determine 
explicit formulas for the spurious exponents in certain cases. Notably, when a 
two-dimensional system with Lyapunov exponents a and (3 is reconstructed in a 
five-dimensional space, we show that the reconstructed system has exponents a, 
(3, 2a, 2(3, and a + (3. 
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Chapter 1 



Introduction 



In the analysis of observed physical systems, it has become standard practice 
to study an auxiliary system reconstructed from a time series of measured data. 
Successful reconstruction of the original system's attractor is the basis of the 
method for Lyapunov exponent calculation proposed by Eckmann and Ruelle 
|l|, U and by Sano and Sawada |§. However, since the reconstructed attractor 
often lies in a larger dimensional space than the original system, the calculations 
produce too many exponents. This leads to an important question: how do we 
distinguish the true Lyapunov exponents of the underlying system from the extra 
"spurious" ones present for the reconstructed system? We answer this question 
for two specific cases: (a) when a one-dimensional system is reconstructed in 
m dimensions and (b) when a two-dimensional system is reconstructed in five 
dimensions. 

The attractor reconstruction process begins by choosing a number m and 
observing the present state p of the underlying system with that number m of 
independent measurements iTi{p), % = 1, . . . ,m. For each point p in the phase 
space, there is an m-dimensional vector n(p) = (vri(p), . . . , n m (p)). This pro- 
duces a measurement function it that associates points in M m with points in 
the phase space of the underlying dynamical system. In practice, the measure- 
ment function 7r often consists of time-delayed versions of a scalar measurement 
(see Takens and others 0, Under certain genericity conditions (see, for ex- 
ample, 0), the original dynamical phase space attractor A will be topologically 
equivalent to its reconstructed image Tt{A) in m-dimensional Euclidean space. 
See Figure 1.1. 

The set of measurement vectors n(p) in the reconstructed attractor tt(A) 
can be studied for geometrical and dynamical properties. Looking to the recon- 
structed attractor for dynamical properties of the original attractor was suggested 
in 1985 by Eckmann and Ruelle, et. al. |l|, 0] and also by Sano and Sawada 0. 
It is usually necessary that m be chosen large enough that there is a one-to-one 
correspondence between points of the original attractor and points of the recon- 
structed attractor ||. This requirement often forces the reconstruction space 
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Figure 1.1: The measurement function tt reconstructs the true attractor A of 
the underlying dynamics (shown in (a)) as the set ir(A) in some M m (shown in 
(b)). When the reconstruction dimension m is large enough, the generic (smooth) 
measurement function n will be a one-to-one correspondence between A and tt(A). 
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to have larger dimension than the original system. In these cases, Lyapunov 
exponent calculations in the reconstruction space produce m real numbers, not 
all of which can be Lyapunov exponents of the original dynamical system. For 
example, if a two-dimensional dynamical system is reconstructed in five dimen- 
sions, current methods compute five "exponents" in the reconstruction space. At 
most two of these can be Lyapunov exponents of the original system; the other 
numbers are "spurious" exponents. 

How, then, do we tell the true exponents from the spurious exponents? Parlitz 
@ proposed that recomputing the exponents using the reversed time series would 
make the true exponents switch sign. Bryant, Brown, and Abarbanel looked at 
the local "thickness" of the data set to identify spurious exponents ||. In addi- 
tion, some authors have proposed removing the extra dimensions by projecting 
the reconstructed dynamics to its tangent plane (see || pp. 336-339] and JT0|, 



pp. 2156-2157]. In the present paper, we study the original Eckmann-Ruelle 
algorithm in order to clarify its output. 

The algorithm presented in || for computing Lyapunov exponents has three 
major steps. First, one reconstructs the attractor in some m-dimensional Eu- 
clidean space of measurements as indicated above. On the reconstructed attrac- 
tor, there is a time-r map F which takes the m-vector P t at time t to the m-vector 
Pt+r at time t + r. This map F represents the reconstructed dynamics, and F(P t ) 
is defined to be P<+ T . In the second step, one computes a local linearization ma- 
trix for F at each point P of the reconstructed attractor by finding the m x m 
matrix M (depending on P) which satisfies as closely as possible 

P t+T - F(P) « M(P t - P) for all t for which P t is close to P. 

We call M = M(P) the Eckmann-Ruelle linearization at P. In the last step 
(which will be discussed in more detail in Chapter 3), one computes the Lyapunov 
exponents of F from a matrix product of these local linearizations. That is, the 
Lyapunov exponents will be the various values achieved by 

lim - In ||M n _iM n _ 2 ■ • • MxM Q u\\ 

n— >oo n 

for various unit vectors v G M m , where Mj = M(Pj) is the best local linearization 
matrix (i.e., the Eckmann-Ruelle linearization) at the point Pj = 7r(pi) in the 
trajectory. 

To understand the output from this algorithm, it is crucial to determine the 
Eckmann-Ruelle linearizations. At first glance, one might think these lineariza- 
tions should produce derivative matrices, DFp. This, however, will not happen 
in general. Suppose that an attractor reconstructed in m-dimensional Euclidean 
space lies within a lower- dimensional surface in M m . The reconstructed dynam- 
ics are well-defined on this surface, but they are not defined off of it, and so the 
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classical mxm derivative matrix DFp will not exist. See Figure 1.1b. Therefore, 
the local linearizations cannot be derivative matrices. 

In this paper, we study the local linearization matrices. For the two cases 
mentioned above, we will show that these matrices have several important prop- 
erties. First and foremost, the linearizations of the reconstructed dynamics are 
surprisingly good in the following sense (see Theorem [Z72] in Chapter 2). For any 
matrix L, if we "linearize" the reconstructed dynamics F about the point P, we 
obtain for some integer k 

F(P t ) - F(P) -L(P t -P)= 0(\\P t - P\\ k ) as P t -> P. (1.1) 

We write g(x) = 0(\x\ h ) as x — > to mean that there exists some constant C 
such that ||p(^)|| < C|a;| fc for all x in some neighborhood of 0. For most choices of 
L, we expect k — 1 in (|1 . 1|) . For the traditional derivative L = DFp, we expect 
k = 2 (when it exists). However, the Eckmann-Ruelle linearization L = M(P) in 
our specific cases has k larger than this. In the case of a one-dimensional system 
reconstructed in m dimensions, we find k = m + 1, and in the case of a two- 
dimensional system reconstructed in 5 dimensions, we find k = 3. The second 
important property of these matrices is that there are natural coordinate systems 
with respect to which the linearization matrices have a special upper-triangular 
matrix representation. This upper triangular form allows us to easily compute 
the Lyapunov exponents (which depend on the diagonal entries when the matri- 
ces are upper triangular). Moreover, for generic measurement functions, these 
exponents will be completely independent of the specific measurement function 
7r used in the reconstruction, depending only on the dynamics of the original 
system. This key property allows us to derive explicit formulas for the Lyapunov 
exponents produced by the Eckmann-Ruelle procedure in the low noise limit. We 
will demonstrate this in Chapter 3. 

Our paper is structured in the following manner. In Chapter 2, we determine 
the linear map that provides the best local linearization to the reconstructed 
dynamics. In Chapter 3, we give the matrix representation alluded to above and 
derive formulas for the Lyapunov exponents produced by the algorithm. The goal 
of Chapter 4 is to prove that numerical determinations of the local linearization 
matrix converge to the Eckmann-Ruelle matrix as the radius of the neighborhood 
shrinks to zero. Chapter 5 presents numerical work illustrating our theoretical 
results. Finally, appendices are included which contain technical lemmas used in 
the text. 
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Chapter 2 



Local Linearizations 



A simple example where f(p) = 2p(mod 2tc) on the interval [0, 2tc] allows us to 
describe the general problem. We reconstruct this interval in M 2 via the measure- 
ment function 7r(p) = (cos(p), sin(p)). The observed dynamics F = vr/7r _1 maps 
points on the unit circle: 7r(p) to ir(f(p)). The following question is central to 
our discussion: for any given point P = 7r(p) on the unit circle 7r([0, 2ir]), which 
2x2 matrix provides the best local linearization of F around P (in the sense of 
fll.l]) )? Since F is not defined off the unit circle, the traditional derivative of F 
does not exist. On the other hand, the linearizations of F are 2x2 matrices while 
the derivatives of the original system / are lxl. It follows that no linearization 
of F can be a derivative matrix. What, then, is the best local linearization of F 
near P (if one even exists)? 

In this chapter, we elucidate the nature of the local linearization matrices 
for the reconstructed dynamics. In their papers jl], Eckmann and Ruelle dis- 
cussed the best local linearization of the reconstructed dynamics. We will define 
the "Eckmann-Ruelle linearization" to be the unique linear map with certain 
properties. In Theorems |2^ and |2l], we will show that our definition provides 
the best local linearization. 

In the simple example above, we specified the measurement function ir ex- 
plicitly. In practice, however, one rarely knows the measurement function that 
arises from the attractor reconstruction process. Time-delay embeddings are 
commonly used, but even then, one cannot know the measurement function com- 
pletely without knowledge of how the scalars in the time-series relate to the states 
of the system. This knowledge is often unavailable in experimental situations. 

For this reason, we focus on measurement functions with generic properties 
(generic in the sense of prevalence) . A property is generic in the sense of preva- 
lence if whenever a function lacks the property, arbitrarily small perturbations 
of that function have the property with probability 1 ||11|| . In the situations that 
interest us, the generic measurement function can be taken to be a smooth diffeo- 
morphism from the underlying phase space into m-dimensional Euclidean space 
M m . The genericity conditions that guarantee the (topological) equivalence of 
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the underlying and reconstructed phase space attractors || typically force the 
reconstruction space to have larger dimension than the underlying attractor. In 
fact, the dimension m of the reconstruction space can be at least twice that of 
the underlying system. 

We adopt the convention that lower-case letters refer to the underlying system 
while upper-case letters refer to the reconstructed system. For example, P = n(p) 
is the point in the reconstructed phase space M m corresponding to the point p in 
the underlying phase space. The convention will also extend to sets in the various 
spaces: u could be a neighborhood of p in the underlying phase space while U 
could be a neighborhood of P in the reconstructed phase space. We represent the 
dynamical flow on the underlying phase space by the time-r map / for some r. 
The measurement function tt maps the underlying phase space into M m for some 
m. In M m , there is an induced time-r flow map F = Trfir -1 that maps each tt(p) 
to 7r(/(p)), and we refer to F as the reconstructed dynamics. Assuming that tt 
is a one-to-one correspondence on the underlying phase space (i.e., that p ^ q 
implies ir(p) ^ 7r(g)), the map F is well-defined on the reconstructed phase space. 

We examine the special case when the underlying map / is one-dimensional 
and the following specific assumptions about /, tt, and a point P = n{p) hold. 

(Al) The underlying dynamical system f maps the unit interval [0, 1] into itself. 
In addition, we assume that f has m+1 continuous derivatives, i.e., f is 

(A2) The measurement function tt maps [0, 1] into M m , and tt is C m+1 . 

(A3) For this p6 [0,1], P = 7r(p) G M m , and the set 7r _1 (P) contains only one 
point, namely p. 

(A4) The first m derivative vectors for tt at p, i.e., 7i^ n '(p) := ^r(p) for n = 
1,2, ... ,m, are linearly independent in M m . 

Properties (A3) and (A4) are generic properties in the space of C m+1 functions 
from [0,1] into M m . By (A4), the derivative vector tt'(p) is nonzero. Together 
with (A3) and the Inverse Function Theorem, this implies 

(A3') There are neighborhoods u C [0, 1] of p and U C 7r([0, 1]) C M m of P such 
that tt(x) G U if and only if x G u, and tt\ u is a diffeomorphism. 

This statement has a useful consequence. Since tt\ u is a diffeomorphism, there is 
a constant C n > 1 such that if Pi, Pi G U with Pj = 7r(pi), Pi G u, i = 1, 2, then 

T^HPi-^H < \pi-p 2 \ <C w ||Pi-P 2 ||. (2.1) 



6 




F(P) 



Figure 2.1: The derivative vectors n^(p) and 7c^(p) (thin arrows) form the 
canonical embedding basis for M 2 at P = n{p). Also, the small tangent vector 
APj = Pi — P maps to its image vector F(Pi) — F(P) near F(P) (thick arrows). 

This inequality will be useful in translating statements back and forth between 
the underlying space and the reconstruction space. 

We begin our analysis of this case by examining the Taylor expansion of n 
about a point P = n{p) in m-space for which (A3) and (A4) hold. Since n is 
one-to-one on the neighborhood u, any point P + AP in ir(u) near P can be 
written in the form 

P + AP = 7r(p + h) = P+J2 + Rem(h) 

n=l n - 

where Rem(h) is the Taylor remainder term (an m- vector here). There is a 
constant Crayior > such that 

\\Rem(h)\\ <C Taylor \h\ m+1 . 

By hypothesis (A4), the vectors n^(p), . . . , 7r^(p) are linearly independent and 
form a basis for ]R m which we call the canonical embedding basis at P. (See 
Figure 2.1.) The little vector AP at P can be written conveniently in this basis: 



AP = 



(h, -h 2 , . . . , ^-h m ) +Rem(h). (2.2) 
V 2 ml J p 



We look at the image F(P) = 7r(/(p)) in essentially the same way. The Taylor 
expansion of F(P + AP) = P(vr(p + h)) = n(f(p + h)) is given by 

m (it o fV^tn) 
F{P + AP) = F(P) + ]T 1 J) {P, h n + Rem f (h). 

n=l U - 
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Here, Rerrif is the Taylor remainder vector associated with 7T o /. Without loss 
of generality, we can choose the constant Cxayior so that we also have 

\\Rem f (h)\\ <C Taylor \h\ m+l . 

At this point, we are ready to define the Eckmann-Ruelle linearization. The 



definition will be justified by the properties stated in Theorem 2.2 



Definition 2.1. Assuming (Al) - (A4), we define the Eckmann-Ruelle lin- 
earization M = M(P) to be the unique linear map from M m to M m satisfying 

M 7r (n) (p) = (TTofY n \p) for each n = l,...,m. (2.3) 

M is well-defined and unique because the set of vectors {tt^(p), . . . , 7i^ m \p)} 
on the left-hand side of fl2.3|) forms a basis for M m by assumption (A4). While it 
may be convenient to think of M as a matrix, no coordinate system has yet been 
specified. We will give a matrix representation for M in the next chapter. Now, 
we prove that M(P) is in fact the best linear approximation to the reconstructed 
dynamics F near P. 

We say that g(x) = 0(\x\ k ) as x — > if there exists a constant C such that 
||g(x)|| < C|x| for all x in some neighborhood of 0. Note that the error term in 
(12. 4) of Theorem can be far smaller than that for the usual Jacobian matrix, 



which would be 0(||AP|| 2 ). 

Theorem 2.2 (Local Linearization, iR 1 ^iR m ). Assume (Al) and (A2). Let 
P = 7T (p) be a point of n([0, 1]) for which (A3) and (A4) hold. If P is in the clo- 
sure of a trajectory of F, P , Pi, • • • e M m , then the Eckmann-Ruelle linearization 
M = M(P) defined by (\2.3\ ) is the unique linear map such that 

F(P + AP) — F{P) — MAP = 0(||AP|| m+1 ) (2.4) 
as \\AP\\ -»• 0, where P + AP G tt ([0,1]). 

Proof. For small AP with P + AP e tt(u), we can write P + AP = n(p + h). 
Note that F(P) = 7r(/(p)) and F(P + AP) = 7r(/(p + h)). First, we prove that 
the Eckmann-Ruelle linearization M = M(P) satisfies (|2.4j ). Using the definition 
of M in equation (p.3|), we cancel terms from the two Taylor series: 



F(P + AP) - F(P) - MAP 

= 7r(/(p + /i)) - 7r(/(p)) - M (7T(p + /i) - 7r(p)) 

/ m w \ / m w \ 

= E^ /)^^)^ + Rem f {h) -ME 7r (n) (p)^7 + ^m(li) 

Vre=l n - J Vn=l n " / 

= Rem f (h) - MRem(h). 
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Since M is a fixed map and \h\ < ||AP|| by ( |2.1| ), we have 

||F(P + AP) -F(P) -MAP|| < || J Rem / (/i)|| + ||MPem(7i)|| 

< c To ^ or (i+iiM|i)|/ i r+ i 

< c Tay/or (i + iiM||)^ +i nAPir +i 



for all ||AP|| sufficiently small, establishing Q2.4|) . 



It remains to show that M is the only linear map that satisfies ( |2.4| ). Suppose 
to the contrary there is another linear map M such that 



P(P + AP) — F(P) — MAP <d\\AP 



I m+1 



as AP — > 0, P + AP G tt(u) 



for some constant C\. Set A := M—M. For small enough AP with P+AP 6 vt(m), 

||AAP|| < ||F(P+AP) - F(P) - MAP 1 1 + ||P(P+ AP) - F{P) - MAP\ 
< C 2 ||AP|| m+1 

where C 2 := Ci + C Ta yior(^ + With respect to the canonical embed- 

ding basis at P, the vector AP has the form in ( |2.2| ), and so, 



A K -h' 



, . . . , 



ml 



< 



AAP\\ + \\ARem(h)\\ 

CTaylor\h\ m+1 



< C 2 \\AP\\ m+1 + 

< C 2 C™ +l \h\ m+1 + \\A\\ C Taylor \h\ m+1 
:= C,\h\ m+1 



Therefore, for h G u sufficiently small 



Mh ,\h*,...,— y h m 

2 ml 



< c 2 3 \h\ 



2m+2 



If we represent the matrix A with respect to the canonical embedding basis at 
P, then the left-hand side is a polynomial in h, call it p(h), with degree at most 
2m. Since P is a limit point of the trajectory in M m , there are infinitely many 
hi G u, P + APj = ir(p + hi) G U, hi — > for which this polynomial satisfies 
< C^\hi\ 2m+2 . By Proposition |A.1| of Appendix A, p(h) must be the zero 
polynomial. This, in turn, forces all the elements of the matrix A to be zero. 
Thus, A = and M = M, proving that M = M(P) is indeed the unique linear 
map satisfying (|2.4|). □ 



To make this method explicit (and for later use), we compute the Eckmann- 
Ruelle linearization for our doubling map example. Let f(p) = 2p(mod 2ir) on 
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[0,27r] and 7r(p) = (cos(p), sin(p)). Let < p < 2tt and set vr(p ) = (x ,y ) 
where x = cos(p ) an d yo = sin(p ). Then, 

* w «=(;ij a,id ^ M= {-2}- 

Also, we have (ir o f)(p) = (cos(2p), sin(2p)) and 

Since Mn^(p ) = (tt o /)«(p ) and M7r( 2 )(p ) = (vr o f)^(p ) by Q, we are 
led to the matrix equation 



M 

which we can solve for M: 



-Vo -%o\ = f-2 (2x y ) -4 (x 2 - y 2 )" 
x -y ) 1 2 (x 2 - yl) -4 (2x y ) 



M(x , y ) 



Axl -Ay, 



2y (2xg + l) 2x (2y 2 + 1) 



In Chapter 5, we will describe numerical experiments for this example showing 
that the linearization matrix calculated by the computer is close to this M(xq, yo). 

We next examine the special case when a two-dimensional underlying system 
/ is reconstructed into five- dimensional space. We make the following specific 
assumptions about /, 71", and a point P = Tr(p) in M 5 . 

(Bl) / maps the unit square [0, 1] x [0, 1] into itself, and f is C 3 . 
(B2) 7r maps [0, 1] x [0, 1] into 1R 5 , and tt is C 3 . 

(B3) For this p G [0, 1] x [0, 1], P = 7r(p) G M 5 , and the set 7r _1 (P) contains 
only one point, namely p. 

(B4) The first and second order partial derivative vectors for tt at p, namely 



^(P) := f^»> %(P) := ^(P) 



^(p) := ^(p), ^(p) := ^(p), T W (P) := ^(p) 

are linearly independent in M 5 . 

Properties (B3) and (B4) are generic properties in the space of C 3 functions from 
[0, 1] x [0, 1] into M 5 . As in the previous case, we obtain a seemingly stronger 
statement as a consequence of (B3), (B4), and the Inverse Function Theorem: 
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(B3') There are neighborhoods u C [0, 1] x [0, 1] of p and U C tt([0, 1] x [0, 1]) of 
P such that tt(x) G U if and only if x G u, and ir\ u is a diffeomorphism. 

This statement implies the obvious analog of ( |2.1| ) for pi G iR 2 and p G iR 5 . 

We begin this case, as before, by examining the Taylor expansions of 7r and 
7T o / about a point P = ir(p) in iR 5 for which (B3) and (B4) hold. With 
h = (hi, h 2 ) G iR 2 , any point P + AP in 7r(w) near P can be written in the form 

P + AP = ir(p + h) 

= P + hiir x (p) + h 2 iTy(p) + -hlir xx (p) 

+ hih 2 TT xy (p) + ~ h^/ffyy (p) + Rem(h,p) 

where Rem(h,p) is the Taylor remainder term, now a vector in iR 5 . Again, there 
is a constant Cxayior > such that \\Rem(h,p)\\ < Cxayior \\h\\ 3 for sufficiently 
small h. The canonical embedding basis (for iR 5 ) at P will consist of the 
five derivative vectors from (B4): 7i x (p), 7l y(p), K xx (p), ir xy (p), and 7i yy (p). With 
respect to this basis, we can write 

AP = [h u h 2 , hi\, hih 2 , + Rem(h,p). (2.5) 

Next, look at the image F(P) = n(f(p)). The Taylor expansion of F(P+ AP) = 
F(iz(p + h)) = n(f(p + h)) is given by 

P(P + AP) = F(P) + h 1 (nof) x (p) + h 2 (nof) y (p) + h 2 1 (7rof) xx (p) 

+ h X h 2 (-K O f) X y(p) + ^hl(lX O f)yy(p) + RCTTl f (/l , Jj) 

Here, Rerrif is the Taylor remainder vector associated with tt o f, and without 
loss of generality, it too satisfies \\Remf(h,p)\\ < CTayior\\h\\ 3 for small h. 

Definition 2.3. Assuming (Bl) - (B4), we define the Eckmann-Ruelle lin- 
earization M = M(P) to be the unique linear map on iR 5 satisfying 

A/r t \ t f\ ( \ Mn xx (p) = (no f) xx (p) 



We will prove in Theorem |2.6| that this linear map M(P) provides the best 
local linearization of F in iR 5 about the base point P. The proof of Theorem 



276| will be virtually identical to that of Theorem |2.2j whose main idea was the 
systematic cancellation of low order terms between the Taylor expansions about 
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a point and its image. Up to second order, there are five distinct terms to be 
eliminated, each requiring a basis element. Thus, we embed into M 5 . Likewise, 
to kill off all third-order terms in a longer expansion, we would need a total of 
nine terms, and so we would embed in M 9 . This is more restrictive than the 
one-dimensional case above since we cannot embed into an arbitrary M m . 

Theorem |2.6| requires a slightly stronger notion of limit point. Since the un- 
derlying phase space has more than one dimension, we should be able to approach 
any underlying base point from more than two directions. This leads to the def- 
inition of an "approach direction." The idea is that there is a line extending out 
from the base point and some subsequence of points from the trajectory approach 
the base point along this line. 

Definition 2.4. Let I be a unit vector in M m . A subset S of M m has the ap- 
proach direction I at the base point P G M m provided there is a sequence 
{QkjkLi from S such that: 

1- Qk — ► P as k — ► oo, and 

2. pQji —>■ I as k — > oo, where £Qk '■= Qk — P ■ 

We call a collection of approach directions at P distinct provided that no two 
are the same and no two are reflections through the origin. Multiple approach 
directions are not required to be linearly independent. 

Figure 2.2 shows three distinct approach directions at a point P on a fractal 
attractor in M 2 . One can find trajectory points converging to P which are on or 
near the intersections of the lines l\, I2, and Z3 with the attractor. This behavior 
was automatic in the one-dimensional case, and points in hyperbolic systems will 
have multiple approach directions. It may even be the case that typical points 
in arbitrary (multi-dimensional) chaotic systems have infinitely-many distinct 
approach directions. 

It is important to be able to relate the (often fractal) geometry of the re- 
constructed attractor back to the geometry of the underlying attractor. One 
reason for this, of course, is that we can only observe the geometry of the recon- 
struction. Another reason is that we need access to facts about the underlying 
system in order to extract information from the observations. These reasons pro- 
vide motivation for the next lemma whose proof is straightforward when it is a 
diffeomorphism of a neighborhood of P. 

Lemma 2.5. Assume (Bl) and (B 2). Let P = n(p) be a point of ir ([0,1] x [0,1]) 
for which (B3) and (B4) hold. If P is in the closure of a trajectory of F, 
Pq, P\, ■ ■ ■ G M 5 that has d distinct approach directions at P , then the underlying 
trajectory po,Pi, ■ • ■ G M 2 , where Pi = 7r(pi), has d distinct approach directions 
at the point p G [0, 1] x [0, 1]. 
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Figure 2.2: The trajectory has three approach directions Li, L 2 , and L3 at the 
base point p. The trajectory is dense in the fractal attractor shown, so data 
points are available arbitrarily close to the intersections of the three lines with 
the attractor. 

We prove the Local Linearization Theorem for the case where we reconstruct 
a two-dimensional underlying system in M 5 . In the statement below, we require 
the reconstructed trajectory in M 5 to have three distinct approach directions at 
the point P. We need three distinct approach directions in order to determine 
the Eckmann-Ruelle linearization. Suppose, for instance, that the available data 
points approach the base point P G iR 5 only along two directions. By Lemma 2.5 , 
the underlying attractor in M 2 will also have two distinct approach directions. 
Suppose that the data lies exactly on the x- and y-axes (in local coordinates in 
M 2 ). We can compute n x (p), ir xx (p), . . . from points of the form (hi, 0). Likewise, 
we can compute n y (p), 7i yy (p), . . . from points of the form (0, h 2 )- However, at 
each of these points, all terms involving mixed partial derivatives vanish from the 
Taylor expansions, making it impossible to compute ir xy (p). Of course, without 
knowledge of ir xy (p) (and where it maps), we cannot uniquely determine a best 
local linearization at P. This is why we require three distinct approach directions. 

Theorem 2.6 (Local Linearization, IB? — ► IB}). Assume (Bl) and (B2). Let 
P = tt(p) be a point of tv([0, 1] x [0, 1]) for which (B3) and (B4) hold. If P is in 
the closure of a trajectory of F, P , Pi, • • • G M 5 that has three distinct approach 
directions at P, then the Eckmann-Ruelle linearization M = M(P) defined by 



6j) is the unique linear map such that 

F(P + AP)-F(P)-MAP = 0(|| AP|| 3 ) (2.7) 
as || AP|| -> 0, where P + AP G 7r([0, 1] x [0, 1]) 
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Proof. We use the same notation as in the proof of Theorem |2.2| , except that we 
write h = (hi, h 2 ) G M 2 . The first thing to do is show that the Eckmann-Ruelle 
linearization M = M(P) satisfies ( J2.7| ). As in the proof of Theorem |2.2| , we use 
the definition of M in equation (|2.6|) to cancel terms from the Taylor expansions. 
Thus, 



\\F(P + AP) - F(P) - MAP\ 



\Remf(h,p) — MRem(h,p)\\ 



< C Taylor {l + \\M\\)\\hf 
for small h. Therefore, by (pT]) , for small AP with P + AP G tt([0, 1] x [0, 1]), 
||F(P + AP) - P(P) - MAP|| < C Taylor (l + ||M||)C 3 || APf , 



establishing ( |2.7| ). 

It remains to show that M is the only matrix that satisfies (|2.7|) . Suppose to 
the contrary there is another matrix M such that 



F(P + AP) - F(P) - MAP <Ci||AP 



as AP -> 0, P + AP G 7r(«) 



for some constant C\. Set A:= M - M. For small AP with P + AP G tt(w), 

||AAP|| < ||P(P + AP) — P(P) — MAP|| + ||P(P + AP) — P(P) — MAP 
< C 2 ||AP|| 3 . 



With respect to the canonical embedding basis at P, the small vectors AP have 
the form shown in (|2.5| ), and so 



.4 [h 1 ,h 2 ,-hl,h 1 h 2 ,^hl 



< 



\AAP\ 



< C 2 1| API 



\\ARem(h,p)\\ 
+ \\A\\ \\Rem(h,p)\\ 



< CzClWhf + WMo 

:= C 3 || /ill 3 



Taylor 



\h\ 



Thus, for all tangent vectors h G u sufficiently small 



A ( h 1 ,h 2 , -hl,hxh 2 , 



<cl 



If we represent the matrix A with respect to the canonical embedding basis at 
P, then the left-hand side is a polynomial p(hi,h 2 ) with degree at most 4. At 
last, we use the assumption that we have three distinct approach directions. It 
follows from Lemma |2.5| , that the underlying trajectory Po,p±, • • • G M 2 also has 
three distinct approach directions. This fact, together with the inequality above, 
gives us precisely the hypotheses needed to apply Proposition |A.3| of Appendix A. 
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Thus, p(hi, h 2 ) must be identically zero, and it follows that all of the coefficients 
of the matrix A must also be zero. Therefore, M = M(P) proving that M(P) is 
indeed the unique matrix satisfying fl2.7|). □ 

This proof extends further to the case when we embed a two-dimensional 
system into higher dimensions. To do this, we first need to embed into an appro- 
priate dimension. In the proof, we cancelled all of the first-order and second-order 
terms in a Taylor expansion. In the general setting, we want to cancel all terms 
of order up to and including order D. There are J2d=i (^t 1 ) = |-^C^ + 3) such 
terms. For each term to be cancelled, we need a basis vector in the canonical 
embedding basis. Therefore, we can embed our two-dimensional system into any 
m-dimensional Euclidean space where m = \D(D + 3) for some D > 0. For 
D = 2 we have m — 5, and for D = 3 we have m — 9. 

We now assume that we are given some D > 1 with m = \D(D + 3) the 
corresponding dimension. We make the following specific assumptions about /, 
7r, and a point P = ir(p) G M m . 

(B5) / maps the unit square [0, 1] x [0, 1] into itself, and f is C D+l . 
(B6) 7T maps [0, 1] x [0, 1] into R m , and n is C D+1 . 

(B7) For this p e [0, 1] x [0,1], P = ir(p) E M m , and the set -k~ x {P) contains 
only one point, namely p. 

(B8) The various partial derivative vectors for Tratpof order at most D, namely 



^XX...Xxip) ■ Q X D (p) 



y > ox ^ ' 

n v(P) := ^(P) 



Txy(p) 



g 2 . . "xx...xx\V) ■ g x D> 

(?) jr ( v ) - _Sf 

d 2 7T / \ l^xx...xy\P) ■— Q x D-i dy \t 

dxdy >* ' ' ' ' 



are linearly independent in M m . 



Properties (B7) and (B8) are generic properties in the space of C D+l functions 
from [0, 1] x [0, 1] into M m . The m vectors in (B8) form the canonical embedding 
basis at P. As before, we define the Eckmann-Ruelle linearization at P to 
be the unique linear map M = M(P) defined by the relations : 

m*m = (w).(p) f/ x f \ = { r f f\ xx ( p \ 

VK J V m ' Mn yy {p) = (7ro/) w (p) 
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Mir xx _ xx (p) = [~ o /).,,,....,,, (/O 

M-Tl xx. ..xy 

(p) = (tt o /) 

xx. ..xy 

. (2.8) 

Af7r w ... w (p) = (7ro/) w .„ w (p) 

Here is the general theorem for two-dimensional underlying dynamics. 

Theorem 2.7 (Local Linearization, iR 2 — > ]R m ). Let m = \D(D + 3) where 
D > 1. ^ssnme f55j and Let P = n(p) be a point of n ([0,1] x [0,1]) 

for which (B7) and (B8) hold. If P is in the closure of a trajectory of F , 
Po, Pi, ■ ■ ■ G M m , and if the trajectory has D + 1 distinct approach directions 
at P, then the Eckmann-Ruelle linearization M = M(P) defined by ( \2.8j ) is the 
unique linear map such that 

F(P + AP)-F(P)-MAP = 0(||AP|| D+1 ) (2.9) 
as \\AP\\ -> 0, where P + AP 6 tt([0, 1] x [0, 1]) 

The proof of this is virtually identical to the proof of Theorem |2.6| and is 
omitted. Note that the polynomial proposition [A.3| is already stated generally 
enough to be used for this proof. 

The Eckmann-Ruelle linearization will also exist in certain other cases, such as 
when reconstructing a three-dimensional underlying system into M 9 . The argu- 
ments above extend in a natural way to certain cases where we have a dynamical 
system / on M n being reconstructed in a larger dimensional space M m . Assume 
that D > 1 and / and 7r each have at least D + 1 continuous derivatives in each 
coordinate. Note that there are J2d=i (^t-i 1 ) = ( D n™) — ^ ^ erms °f order at 
most D in the Taylor series of / and ix. Then, when m = ( D ^ n ) — 1, we can 
construct the canonical embedding basis as in assumption (B8) and define the 
Eckmann-Ruelle linearization M(P) as in ( |2.8| ). Because ( p.8|) guarantees that 
the Taylor series will collapse nicely, it is easy to see that M(P) satisfies (|2.9|). 
We will also need some mild geometric condition similar to those given in the 
previous theorems to guarantee that M(P) is the only matrix satisfying (|2.9| ). 
We will not pursue these ideas further at this time. 
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Chapter 3 



Lyapunov Exponent Formulas 



Once the local linearization matrices have been computed at each point of the 
reconstructed trajectory, we must extract Lyapunov exponents from them. Mim- 



icking the standard definition of Lyapunov exponents (see for example, fll|, pp. 
31]), we define the Eckmann-Ruelle-Lyapunov (ERL) exponents of the re- 
constructed trajectory P , Pi, ... in M m to be the values obtained by the limit 

h ER (Po, v) ■= lim - In \\M n ^M n „ 2 ■ ■ ■ M 1 M u\\ (3.1) 

for unit vectors v G M m , where M{ = M(Pj) is the best local linearization (i.e., 
Eckmann-Ruelle linearization) at the point Pi = n(pi) in the trajectory. We need 
a separate definition here because the standard definition of Lyapunov exponent 
uses the Jacobian derivative, which need not exist along the trajectory. In this 
chapter, we show that the matrix product in ( |3.1| ) can be written as an upper 
triangular matrix. A straightforward calculation will then produce a formula for 
the limiting values of ( |3.1| ). 

In practice, the limit in (3.1) can be computed using the treppen-iteration 
algorithm described in [|T], Q and elsewhere. This method uses QR matrix de- 
composition to convert the matrix product M n _i • • • M into a product of upper 
triangular matrices R n -i ■ ■ ■ Ro- The diagonal elements of the latter upper tri- 
angular matrix are the products of the corresponding diagonal elements of the 
Ri. Then, we can read the Lyapunov exponents right from the diagonals of the 
intermediate matrices Ri. 

n-l 



A,- = lim — In \ (Ri 

J n— >oo <n ^ — 4 
H i=0 



J J 



Theorem p. 2| justifies the ability to read the exponents directly from the diagonal 
in this way. A proof of the theorem will be given in Appendix B. 

Definition 3.1. A sequence of real numbers {r n } has (geometric) growth rate 

7 provided 



In \ r 



lim — — — = 7 



-•oo 



n 
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Theorem 3.2. For k = 1,2..., let be m x m upper triangular matrices, 
and define S n = A n ---A 1 . Assume the magnitudes of the entries of A k are 
bounded independent of k, and that the diagonal entries of S n have growth rates 
71, . . . , 7 m as n — > oo. Then there exist vectors v%, . . . , v m G M m such that for 
each % — 1, . . . , m, \\S n Vi\\ has growth rate 7*. 



The next theorem and its two-dimensional analog, Theorem |3.6| , are the main 
results of this paper. In both theorems, most of the hypotheses are used to guar- 
antee that the Eckmann-Ruelle linearization exists at each point of the trajectory. 

Theorem 3.3 (Lyapunov Exponent Formula, JR 1 — > ]R m ). Assume (Al), 
(A2) and that the trajectory of f in [0, 1], Po,pi, ■ ■ ■ , has Lyapunov exponent X. 
If each point of the reconstructed trajectory, Pi = 7r(pj) G M m , satisfies (A3) and 
(A4) and is a limit point of the trajectory, then the reconstructed trajectory has 
Eckmann-Ruelle-Lyapunov exponents A, 2 A, 3 A, . . . , mX. 



In order to prove this, a few lemmas will be useful. Lemma |3.4| is a technical 
lemma leading to the matrix representation for M(P) given in Lemma |3.5| . With 
this matrix representation, we will prove Theorem |3.3| . 



Lemma 3.4. Let J be an interval, and let k,d be positive integers. If f : J —> J 
is a C k+1 function, then for 1 < j < i < k there are differentiate scalar functions 
dijit), such that for any C k+1 function (ft : J — > M d , <f>(t) = (<fii(t), . . . , 4>d(t)), we 
have 

(<P o ff it) = £ a l3 (t)^\f(t)), fori = l,...,k (3.2) 
j'=i 

where (f)^\t) = (^^r(^), • • • , ^^(^)) represents the j-th derivative vector of <f>(t). 

Note that the scalar functions aij(t) depend only on / and not on the other 
function <f>. As will be seen in the proof, these functions are the results of 

collecting terms involving the derivatives of <fi. 



Proof of Lemma |3.4| . Since vector-valued functions of a real-variable are 
naturally differentiated componentwise, it is enough to prove the lemma for any 
one component. Thus, we may assume d — 1 and <fi : J — > M. 

The proof is by induction. For % = 1, we apply the Chain Rule: 

(t)=f(t)0«(/(t)) 
and note that an(t) := f'(t) is C k . For i = 2, we apply the chain rule again: 
(0 o /) (2) (t) = /'(t) 2 (2) (/(*)) + f" (W {1) (f(t)). 
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Note that a 21 (t) := f"(t) and a 22 {t) := {f'{t)) 2 are both C k ~ l . 

For induction, suppose that ( |3.2|) holds for i where the are C fe+1 ~ l , 

depending only on /, not on 0. We differentiate ( |3.2| ) and collect together terms 
involving the derivatives of <fi. Specifically, if we define the functions au+i)j(t) 
using 

:= a' a (t) 

a (i+i)j(t) a 'ij(t) + a iti-i)( t )f'( t )> j = 2, . . . , z— 1 

:= i7'(*) i " 1 /"(0 + a<(i-i)(*)/'(*) (3-3) 

a(t+i)(t+i)(o ;= f(ty + 

then we obtain from differentiating ( |3.2|) : 

(0 ° ff +1) (t) = j t [f'm {i) {f{t)) + 

= 5><m)iW W (/(*))- 
j'=i 

Note that the functions a^ +1 )j(t) are (7 fc +i-(M-i) anc l that they depend only on 
the function /, not on the other function <fi. This completes the proof. □ 

Lemma 3.5. Assume (Al ) and (A2). Let P = n(p) G M m and let Q = F(P) be 
points satisfying (A3) and (A4) in the closure of a trajectory of F. Let a and (5 
represent the canonical embedding bases at P and Q, respectively. With respect to 
these bases, the Eckmann-Ruelle linearization [M(P)]@ is upper triangular with 
diagonal elements M(P)jj = f'(p) J ■ 

Proof. Recall that the canonical embedding basis at a point P = ir(p) in 
lR m consists of the vectors n^{p), . . . , Tr (m \p) in R m . Set M = M(P) and 
q = f(p). By the definition of the Eckmann-Ruelle linearization (|2.3|) , we have 
Mtt^(p) = (it o fp l \p). We must write the vectors (it o f)^(p) in terms of the 
canonical embedding basis at Q = ir(q). Applying Lemma [3.4| with k = d = m 
and (ft = 7r, we have 

M^\p) = (tt o /)« (p) = /'(p)V%) + £ <k> (p)* U) {q) for i = 1, . . . , m 

3=1 
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because f(p) = q. This provides a representation for M with respect to the 
canonical embedding bases a and (3: 





(f(p) 





02l(p) 

f(p) 2 




031 (?) ■■ 

032 fa) •■ 

/'(p) 3 


■ oh(p)^ 

■ afe2(p) 

■ afe3(p) 




{ 








• /'(p) m y 



It is important to note that this matrix form for the Eckmann-Ruelle lin- 
earization is completely independent of the embedding it. Information about it is 
used to form the coordinate system with respect to which we view the dynamics; 
however, once inside that coordinate system, we only see the action of /. 

Proof of Theorem |3.3| . For each point Pj e M m , the Eckmann-Ruelle lin- 
earization Mi = M(Pj) exists. To compute the ERL exponents, we must evaluate 
the limit (|3.1| ). Let j3i denote the canonical embedding basis at Pj. Notice that 
the canonical embedding bases for the upper triangular matrix representation fit 
together perfectly: 

It follows that the matrix product M rt _ 1 M n _ 2 ■ • • MiM Q is upper triangular when 
expressed with respect to the canonical embedding bases (3 and (3 n . That is, the 
matrix product [M n _ 1 M n _ 2 ■ ■ -MiMq]^ can be written 



/n-l 

n f(Pi 

3=0 





'12 




>2?, 




'Ik 



>2k 



'n-l 



V 



n f(pj) 



J 



where the fry are numbers which depend solely on the underlying dynamical 
system / and the first n points of the trajectory of po in [0, 1]. 

We apply Theorem |3.2| to complete the proof. Note that the elements of each 
Mi (with respect to the appropriate coordinate systems) are combinations of the 
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first m derivatives of /, each of which is continuous. These elements will be 
bounded independent of i on any compact set containing the entire trajectory. 
By hypothesis, the Lyapunov exponent A of / is given by 



A = lim — In 

n— >oo yi 



n— 1 

n fipi. 

3=0 



and it is clear that the diagonal entries of the matrix product M n _! • • • M 



have growth rates A, 2A, . . . , mX. Theorem |372] now provides the vectors which 
grow at these characteristic rates, and we conclude that the ERL exponents are 
A, 2A, ...,mX. □ 

With this theorem, we expect that whenever the underlying system is one- 
dimensional with Lyapunov exponent A, we will compute exponents A, 2A, . . . , mX 
in the absence of noise. Note that if the system has a positive Lyapunov exponent 
(and hence is chaotic), then the true exponent will be the smallest of the computed 
numbers! 

Next, we consider the case where the underlying system is two-dimensional 
and reconstructed in M 5 . The basic arguments are similar though the situa- 
tion is more involved, as we shall see. The formula for the computed Lyapunov 
exponents in this case is given in the next theorem. 

Theorem 3.6 (Lyapunov Exponent Formula, M 2 — > M 5 ). Assume (Bl), 
(B2), and that the trajectory of f in [0, 1] x [0, 1], po,pi, ■ ■ ■ , has Lyapunov ex- 
ponents X and \x. Assume each point of the reconstructed trajectory, Pi = n{pi) 
in M 5 , satisfies (B3) and (B4), and the set {Pq> Pi, P21 ■ ■ ■ } M 5 has at least 
three distinct approach directions at each Pi. Then, the reconstructed trajectory 
has Eckmann-Ruelle-Lyapunov exponents X, fi, 2X, X + fi, and 2fi. 



Proof. Recall that the canonical embedding basis given in assumption (B4) 
from Chapter 2 consists of the first and second order partial derivative vectors of 
7r, namely 7T x (p), n y (p), n xx (p), ir X y(p), and TT yy (p). In fact, the construction of 
the Eckmann-Ruelle linearization from the Taylor series of / and ir can be carried 
out with respect to any orthonormal set of coordinates. The uniqueness part of 
the Local Linearization Theorem |2.6| ensures that the resulting linear map will be 
the same. Thus, we may introduce convenient local coordinate systems at each 
point Pi G M 2 of the trajectory of p . 

Without loss of generality, we assume A > \i. We begin by choosing a unit 
Lyapunov vector vq G M 2 corresponding to the exponent fi at the initial point 
Po G M 2 of the underlying trajectory. That is, we choose a unit vector vq such 
that 

1 



lim — In 

n^oo jy 



Dfpn^Dfpn-2 ■ ■ ■ Df Pl Df Po vo 
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By Oseledec's Multiplicative Ergodic Theorem |]J, almost every other vec- 
tor in M 2 has growth rate A. Choose a unit vector wq perpendicular to vq. 
This gives us an orthonormal basis {vq,Wq} for M 2 based at p . In general, 
given the basis {v n ,w n } at p n G M 2 , we construct the basis at p n +i by setting 
v n+ \ := Df Pn v n /\\Df Pn v n \\ and choosing w n+ i to be the unit vector perpendicular 
to v n+ i that satisfies (w n+ \, Df Pn w n ) > 0. Thus, at each point p n , we have an or- 
thonormal basis for M 2 . For each n, we write points p' G M 2 near p n as p' = (x, y) 
provided p' = p n + xt> n + yw n . Now, we can describe the underlying dynamics 
/ near p n in terms of these bases at p n and p n +i by f(x,y) = (g(x,y),h(x,y)). 
Though we will not explicitly show it, this representation for / depends on the 
base point p n , and the reader should keep in mind that the component functions 
g and h may look very different as we vary p n . Note that h x {p n ) = because 
Df Pn v n = \\D f Pn v n \\v n+ i has no y-component at p n +i- Thus, we can write the 
Jacobian derivative of / at p n as: 



Df p 




Following the ideas in the previous case, a calculation using the chain rule pro- 
duces the next set of equations as the analog of (|3.2|), where we write q n = f(p n ) 
and the partial derivatives of g and h are evaluated at p n : 



° f)x(Pn) = 9xK x (q n ) 

(n°f) y (Pn) = g y ^x{q n ) + hy-K y (q n ) 
(tt o f) xx (p 

n) gxx^x (?n) + h xx n y (q n ) + g 2 x 7i xx {q n ) 

(7T O f) X y(p n ) = g X y1T x 

O f)yy{p n ) = g yy TC X (q n ) + hyyir y (q n ) + ^ 7T XX ( <?„ ) + kyTT X y fa) 



hl^yyfa) 



For each n, let /3 n denote the canonical embedding basis for M 5 at P n = 7r(j9„). 
Representing the Eckmann-Ruelle linearization with respect to (3 n and (3 n+ i as 
we did in Lemma [3.5| , we obtain this matrix representation for M n = M(P n ): 



fax 
o 
o 





9y 
h y 









gxx 

hxx 

9l 





9xy 
foxy 
9x9y 

g x h y 





9yy 

h. 



\ 



2 



9y 

2g y h y 
h 2 



All of the terms above the diagonal come from combinations of first and second 
derivatives of /, and, because / is C 3 by hypothesis, they will be bounded inde- 
pendent of p n . Note that the upper left 2x2 block is the Jacobian Df Pn of the 
underlying dynamics and that the lower 3x3 block contains only combinations 
of terms from Df Pn . 
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As in the proof of Theorem |3.3| , the product [M„_i • • ■ M ]^ is upper trian- 
gular when written with respect to the canonical embedding bases: 

/n-l \ 

II 9x(Pj) &12 &13 & 14 615 

3=0 

n-l 

n ^(Pj) & 23 &24 &25 

fjI^(Pi)j b M 6 35 

n-l 

II 9x(Pj)h y (pj) 645 

i=o 

(n-l 
n 
i=o " J J 

where the 6jj are numbers which depend solely on the underlying dynamical 
system / and the first n points of the trajectory of p in [0, 1] x [0, 1]. Now, we 
need to compute the growth rates of the diagonal terms. To do this, we compute 
the Lyapunov exponents A and fi of the underlying dynamical system in terms 
of the components of the Jacobian matrices. Since Df Pn v n = \\Df Pn v n \\v n+ i, we 
have g x {p n ) = \\Df Pn v n \\. It follows that 

lim - In ( TT g x {p n )) = lim -In ( TT \\Df Pt Vi\\) 

1, „„, „ (3-4) 



£™o- ln \\ D fpni--- D f P o V o\ 



because vq was chosen to be a Lyapunov vector for \i. Recall from JI], pp. 632] 
that the growth rate of areas is given by the sum of the Lyapunov exponents. 
Thus, 



+ \ — Jim - In |det (Df n po , 



n— »oo fi 

/n-l 



1 ( n ~ l 
lim - In TT \g x {pn)hy{p r , 
n Vl= o 

/i + n lil ^~ ln (S \ h v(Pn)\ 
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and it follows that 

A= ii^^ n (niwi)- (3.5) 

Finally, we see that the diagonal terms have growth rates of A, fi, 2 A, A + /x, 
and 2fi. It follows from Theorem [3.2| that the reconstructed trajectory, P , Pi, . . . 
in M 5 has these ERL exponents, completing the proof of Theorem |3.6| . □ 
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Chapter 4 



Convergence Theory 



In Chapter 2, we showed that, at least in theory, the local linearization matrices 
are not derivatives. The formulas derived in Chapter 3 for the Lyapunov expo- 
nents are meaningful only if the numerically determined linearization matrices 
are close to the Eckmann-Ruelle linearizations M(P). Because M(P) is guaran- 
teed to be the best linearization in the limit as the neighborhood radius shrinks 



to zero (see Theorems |2.2| and |2.6|) , it is possible that for a particular radius, a 
different matrix will do "better" than the Eckmann-Ruelle linearization, though 
we hope the "better" matrix will still be close to M(P). In this chapter, we prove 
that, over small neighborhoods around the base point P, the best linearization 
matrix will in fact be close to M(P). 

To proceed, we must specify what we mean by the "best" linearization over a 
small neighborhood. That is, we need a means of measuring the error involved in 
the linearization process. The following definition serves this purpose by finding 
the worst-case error committed by a matrix L used as the local linearization. 

Definition 4.1. Given F : lR m — > lR m , an trajectory of F, P ,P ir -- e M m , 

and P G M m in the closure of the trajectory. For an m x m matrix L and e > 0, 
define 

W{L, P, e) := sup ||P i+1 - F(P) - L(P t - P) \\ 

\\Pi-P\\<e 

where the supremum is taken over those values of i for which \\Pi — P\\ < e. For 
matrices L\ and L-2, we say that L\ is a better linearization than L 2 over the 
e-ball about P provided W{L\, P, e) < W(L 2 , P, e). 

In practice, one needs to avoid false-nearest-neighbors by taking the supre- 
mum over % where both \\Pi — P\\ < e and ||Pj+i — F(P)|| < e. However, this is 
not necessary for the theory because 7r is one-to-one in a neighborhood of p (by 
either (A3') or (B3')) and we will let e — > 0. Hence, we can assume e is small 
enough that the false-nearest-neighbor problem never arises. 

We point out that the matrix minimizing W(-,P,e) is not necessarily the best 
least squares fit but rather the best "minimax" fit. The least squares problem 
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seems harder to formulate theoretically. At the end of this chapter, we present 
preliminary results for a least squares theory that incorporates infinite data. 

In the statements of Theorems ^72| and [O] below, we assume that we have a 
sequence of "best" linearizations M k over a shrinking set of e^-balls. This is the 
interpretation we give to the condition W(Mk, P, e k ) < W(M(P), P, e k ). Thus, 
each matrix M k does better than M(P) over the e^-ball around P though not 
necessarily over any other ball around P. 

As before, we begin with the one-dimensional case where / : [0, 1] — > [0, 1] 
and it : [0, 1] — > iR m . We say the reconstructed trajectory is dense in the 
curve at the point P G 7r([0, 1]) provided there is a neighborhood U of P such 
that the trajectory {P : i > 0} is dense in the curve 7r([0, 1]) fl £/. 

Theorem 4.2 (Convergence Theorem, M 1 — ► _K m ). Assume (Al) and (A2). 
Let P G 7r([0, 1]) fre a pomt satisfying (A3) and (A4) where the reconstructed 
trajectory is dense in the curve. Let {tk}kLi be a decreasing sequence of positive 
numbers, e k — > 0. Let {M^}^ be a sequence of matrices for which we have 
W(M k , P, e k ) < W(M(P), P, e fc ) for all k. Then M k -> M(P) as fc -> oo. 

The key point for the proof is to control the size of \\M k — M(P)\\ as e k — > 0. 
This control will be achieved via the next proposition. Recall that a convex set 
contains every line segment connecting any two of its points. That is, S is convex 
provided that whenever u,v G S and t G [0, 1] we have tu + (1 — t)v G S. The 
convex hull of a set S is the smallest convex set containing S. 

Proposition 4.3. Let S be a subset of M m , and let Hull(S) denote its convex 
hull. Assume that Hull(S) contains a closed m- dimensional ball of radius r. If 
A is an m x m matrix such that \\Av\\ < B for all v G S, then \\A\\ < 

Proof. Note that \\Av\\ < B actually holds for all v G Hull(S). Let rj be a unit 
vector such that \\Ar]\\ = \\A\\. Since Hull(S) contains the closed ball B(c,r) for 
some c G Hull(S), there exists z G Hull(S) on the surface of this ball for which 
the vector z — c points in the direction of rj. Therefore, z — c = rt], and we have 

r \\A\\ = r \\At]\\ = \\A(z - c)|| < ||As|| + ||Ac|| < 2B. 

□ 

Here is how we will use the proposition. We let A k be the matrix M k — M(P) 
and S k the set of small tangent vectors from the efc-ball about P. The bound 
on ||Au|| will come from the condition W(M k ,P,e k ) < W(M(P),P,e k ). We will 
then be able to conclude a bound on \\M k — M(P) || in terms of e k . In order to 
proceed, then, we must examine the convex hull of the tangent vectors. 
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Figure 4.1: The solid curve outlines the convex hull of the curve (h,\h 2 ) for 
h G [—1, 1]. No small perturbation V(h) (dashed curve) can significantly change 
the interior of this convex hull. Thus, the convex hull of V(h) will always contain 
a ball of radius C ra d . 



Recall that when written in the canonical embedding basis at the base point 
P, little tangent vectors AP, := Pi — P have the form shown in (^2). This 
suggests that we consider curves in M m of the form 



V(h) = [h,-h 2 , 



—h m + Pert(h) 
ml 



where Pert(h) is a continuous vector-valued function with Pert(O) 
Hulljy) denote the convex hull of the set {V(h) : h 6 [— e, e]}. 

We need some additional notation. For i = 0,1, ... ,m, define x$ 
Vand m denote the m x m Vandermonde matrix: 



(4.1) 
0. Let 
■i-. Let 



Vandr, 



/I 
1 



X'2 



\1 X r , 



x\ 
.}, ■-) 



m— 1 
X% 



X 



m—1 



Since determinants vary continuously with the elements of the matrix, there is 
eo > such that if we perturb the elements of Vand m by no more than eo then the 
determinant of the perturbed matrix, Vand m + E, is still close to the determinant 
of the unperturbed matrix Vand m : 

if < e o f° r au i an d j, then det (Vand m + E) > — det (Vandm) > 0. 
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Lemma 4.4. There exists a constant C m <z > such that for any curve V(h) of 
the form U-A ) satisfying max_i</,<i ||Pert(/i)|| < ^f^j, the convex hull HulliiV) 
contains a ball of radius C ra d- 



Proof. The m + 1 points V(xi) in M m , i = 0,1, ... ,m, form an m-simplex 
entirely contained (along with its interior) inside Hulli(V). Note that V(xq) is 
the zero vector. Treating the V(xi) as row-vectors and writing £y for the jth 
component of Pert(xi), the volume of this simplex is given by [[i3j pp. 44]: 



det 



V(x 1 ) 




Xl + Cu 


V{x 2 ) 




^2 + 6l 


= — det 






ml 




V(x m ) 







2 X 1 + 6-2 
2^2 + £22 



m \Xl 4" ^lm 



2 3 'r?i 4" £' 



)/i2 



1 ^.m 
m! m 



ml 



ml m' 



n 

\A:=1 



det 



1 + ^^11 
1 + ^21 



xi + ^12 

^2 + —62 



m—l 



> 



m" 




1 + 

det (Vcmc/ m ) > 



ml 



X 2 



m—l 1 mi 



Irn 



+ 



because the perturbations of the elements of Vand m satisfy: 



ml 

< — < m m! II Pert(xi) || < e . 

X\ 



Since the vectors V^Xj) cannot vary too far from their "unperturbed positions" 



at the points [x^ ^xf , . . . , ^jX^J, and since the m-simplex they form has an 
absolute lower bound on its volume, there must be some radius C ra( i > such 
that this m-simplex always contains a ball of that radius, no matter what the 
perturbation function Pert(h). □ 

Note that the vectors V(xi) and the Vandermonde matrix used in this proof 
produce a simplex with very small volume. Indeed, other choices for the Xi could 
produce larger volumes and hence larger "guaranteed" radii C ra( i- Fortunately, 
we do not need the radius to be large, only that such a guaranteed radius does 
in fact exist. 

Proposition 4.5. Let Cp ert > and < e < min (l 



' m m\Cp e 



For any curve 

V{h) of the form g^) with \\Pert(h)\\ < C Pert \h\ m + l defined on \-e,e], Hull € (V) 
contains a ball of radius C ra d£ m , where C ra d is independent of e and V(h). 
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Sketch of Proof. This is a consequence of changing coordinate systems 
between the standard basis {ex, e 2 , . . . , e m } and the basis {eei, e 2 e 2 , . . . , e m e m }. 
Then, as h ranges over [— e, e], - ranges over [—1, 1]. Together, the conditions on e 
and Pert(h) imply that the perturbation term is still small enough to use Lemma 
[4.4| . This guarantees the existence of a ball of radius C ra d contained within the 
convex hull of V over [—1, 1]. Naturally, this ball has axes of length C ra d in 
each of the principal directions of the second basis. When we convert coordinates 
back to the standard basis, the ball becomes an ellipsoid with axes of length 
Crade 1 for i — 1,2, ...,m. In particular, the ellipsoid contains a ball of radius 

C ra d£ m - □ 



Proof of Theorem |4.2| . The hypotheses guarantee that the Eckmann-Ruelle 
linearization M(P) exists and is the unique best linearization as e k — ► 0. We will 
show that there is a constant C (independent of k) such that ||M fc — M(P)\\ < Ce k 
for large k. 

By Theorem [2.2| , there are e' > and C er > such that for every < e < e' 
we have W(M(P),P,e) < C er e m+1 . Without loss of generality, e% < min{e', 1}. 
Then, 



W(M k ,P,e k )<W(M(P),P,e k ) 



771+ 1 



for k = 1, 2, 



(4.2) 



First, we construct a suitable set S k . With respect to the canonical embed- 
ding basis at P, AP; := P t - P = (h u \h\, ^hf) + Rem(hi). Without loss 
of generality we may assume that Rem(h) is defined on [— ei, ei]. Define the func- 
tion V(h) on [— ex,ei] by equation ( |4.1| ) with Pert (hi) = Rem(hi), and note that 
APj = V(hi). Since the reconstructed trajectory is dense in the curve at P, we 
may assume that the points Pj in B(P, ei) trace out the curve 7r([0, 1]) fl B(P, ex). 
Because ir is a diffeomorphism, these points Pi G M m pull back to scalars pi G M 
that are dense in some subinterval of [0, 1]. By fl2.1|) , this interval contains the 



subinterval 
interval Ii 
P % G B(P,e k )nir 



a 

ei 

C 7r 



It follows that the /i,- 



Pi 



p are dense in the 



ei 

Ctt 



By similar reasoning, for each /c, the preimages of 



1]) are dense in the interval p — + jf- 



and the corre- 



sponding hi are dense in I k 



6fe 

^7T 



Define the set 5*^ := {V(h) : /i G I k } in 



iR m . The convex hull of S k is Hull^k_(V), and by Proposition [4.5| , it contains a 
ball of radius ^rejT whenever < min ^1, 



m m!C TaI/ior . 

For each k, let := M k — M(P) be the m x m error matrix. Then, for each 
i with ||APj|| < e k , we have by 



\A k AP t \\ < \\P i+ i - F(P) - M(P)APi\\ + \\P i+1 
<W(M(P),P,e k ) + W(M k ,P,e k ) 



F(P)-M k APi\ 



< 2C er e k n+l 
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Because APj = V(hi), the hi are dense in 1^, and V(h) is continuous on this 
interval, this inequality extends to all h G Ik- Therefore, 



\\A k v\\ < 2C er e™ +i for all v e S, 



k ■ 



and, by Proposition [O 



HA IK ACerC * c 

\\ A k\\ — e k- 

This completes the proof. □ 
Before we proved Theorem [4.2| in the one- dimensional case, we added an as- 



sumption about the distribution of the orbit within the reconstructed attractor 
that was not needed in previous chapters. In particular, we assumed that, in a 
neighborhood of the base point, the orbit was dense along the curved segment 
of the reconstructed attractor. Unfortunately, a density assumption is much too 
strong for the case of two-dimensional underlying dynamics because our recon- 
structed attractor could be fractal in nature. 

It turns out that we need only add a third condition to the definition of 
"approach direction." This condition ensures that the distances between the base 
point and points approaching along the direction vector shrink exponentially fast. 



Definition 4.6. Let I be a unit vector in M m . A subset S of M m has the fractal 

i oo 

lfc=l 



approach direction I at the base point P G M m if there is a sequence {Qk}^ 



from S such that: 

1- Qk —> P (is k — > oo ; 

2. pQji I k — > oo ; where £Qk '■= Qk — P, and 

3. there are constants < C~ < C + < 1 such that for all k 

C~ < < C + . (4.3) 

" II^Qfcll ~ 

A collection of fractal approach directions at P is distinct provided that no two 
are the same and no two are reflections through the origin. Multiple fractal ap- 
proach directions are not required to be linearly independent. 

As noted in Chapter 2 (see Figure 2.2), this kind of behavior occurs for points 
in hyperbolic systems and may occur in all chaotic systems with two or more 
dimensions. 



We need a lemma like |2J| that translates fractal approach directions in the 
reconstruction space down to fractal approach directions in the underlying phase 
space. 
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Lemma 4.7. Assume (Bl) and (B2). Let P = n(p) be a point ofn([Q, 1] x [0, 1]) 
for which (B3) and (B4) hold. If P is in the closure of a trajectory of F, 
Po, Pi, ■ ■ • G M 5 that has d distinct fractal approach directions at P, then the 
underlying trajectory Po,Pi, •• • G M 2 , where Pi = ir(pi), has d distinct fractal 
approach directions at the point p G [0, 1] x [0, 1]. 



Proof. By Lemma |275|, we get d distinct approach directions in M 2 . Thus, 
we need only verify the third condition for each one. Let / be an approach 
direction in M 2 coming from a fractal approach direction in M 5 . Let Qk = ^{qk), 
k = 1, 2, . . . , be the points in M 5 forming the subsequence that converges to this 
fractal approach direction. By the previous comment, q k — > p and iffiE^] ~~ * I- We 
can use the constants C~, C + from condition (3) together with (|2.1|) to translate 
(3) to: 

c_ ||g fc+1 -p|| < c+c2 
c 2 hk-p\\ 

If C + C 2 < 1, we are done. Otherwise, we have to come up with new bounds. 
We form a new subsequence from {g^jfeLi as follows. Start with qi = q%. Then, 
having picked qk = qi for some i, we must pick q~k+i- Look back to the original 
sequence {qk} and set qt+i to be the first term after q iy say qi + j with j > 0, for 
which \\qi + j — p\\ < C + ||gj — Note that ||gj +3 -_i — p\\ > C + ||gj — p\\- Then, 

c + > Mk+i -p\\ = H+j-pW H+j-i -p\\ > 9_ c + 

Il?i+J-1-P|| hi~P\\ ~ C l 

These bounds now satisfy (3) of the definition. Of course, since {%} is a subse- 
quence of {qk}, we have qk — > p and ii^^ii - ► I as well. □ 

The Convergence Theorem below applies only for the case of a reconstruction 
in M 5 , but the extensions to appropriate higher- dimensional cases should be clear. 

Theorem 4.8 (Convergence Theorem, M 2 — > IB}). Assume (Bl) and (B2). 
Let P = 7r(p) be a point of tt([0, 1] x [0, 1]) for which (B3) and (B4) hold. Assume 
that P is in the closure of a trajectory of F, Pq,Pi,-- - G M 5 that has three 
distinct fractal approach directions at P. Let {e^}^ be a decreasing sequence 
of positive numbers, e& — ► 0. Let {M^^i be a sequence of matrices such that 
W(M k , P, e k ) < W(M(P), P, e fc ) for all k. Then M k -> M(P) as k -> oo. 

This proof will be similar to that given in the one-dimensional case over the 



course of Lemmas |4.3| , |4.4| , and |4.5| . The present situation will be far more techni- 
cally involved, however, because the geometry is no longer simple. In the previous 
case, we had only a single approach direction and it was reasonable to assume 
the trajectory was dense along the line. These properties are not generic when 
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Figure 4.2: For this configuration of points z%, . . . , z§ in M 2 , the convex hull in 
M 5 of the set {tt(zx), . . . , vr(z 5 )} is guaranteed to contain a ball of radius C ra d- 



the underlying dynamics have two or more dimensions. (See Figure 2.2.) In the 
present case, the local distribution of the trajectory can be sparse, approaching 
only along three fractal approach directions. Recall from the discussion in Chap- 
ter 2 that three approach directions are necessary to ensure that we can uniquely 
determine the Eckmann-Ruelle linearization. 



In Lemmas [O] and |4.5| , we singled out a configuration of points in the under- 
lying space (easily chosen thanks to the density assumption) whose reconstructed 
images had a convex hull that could be guaranteed to contain a small ball. We 
will do the same thing here, though we have to be more careful. Our points 
are no longer nicely constrained to a line (they converge along three approach 
directions in the plane), and we cannot choose any points we want (there is no 
density assumption). We will choose five points Zi situated near our three fractal 
approach directions U at P as in Figure 4.2. We will use the spacing between 
points z% that (|4.3|) guarantees in the same way that we used the guaranteed 
spacing between x.i = — and = — in the one- dimensional case. 

To be more specific, for i = 1,2,3, let U = (cos(aj), sin(aj)) be the fractal 
approach directions. Without loss of generality, we can label the U to satisfy 

< ai < 2tt for i = 1, 2, 3 

< a 2 - ai < it (4.4) 
< a 3 — a 2 < 7r 



Each li will have associated constants C i and from (|4.3|) , but we can easily 
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find constants < C~ < C + < 1 that work for all three Zj simultaneously. We 
want to consider quintuples of points z$ = (rj cos(^), r» sin(0j)) G iR 2 , 1 < i < 5, 
that satisfy these constraints for some rj > 0: 



C-<r,<l i = l,3,4 If 1 "" 1 !-^ Z = 1 ' 2 

r ' ' — a 3 | < 77 * = 4,5 (4.5) 

We will also consider functions V : iR 2 — > M 5 of the form 

V{x, y) = (x, y, ^x 2 , xy, ^y 2 ^j + Pert(x, y) (4.6) 

where Pert(x, y) is a continuous vector- valued function. In the following lemma, 
the function V plays the role of our measurement function n, written in canonical 
embedding coordinates. 

Lemma 4.9. Given constants < C~ < C + < 1 and ai,a2,as satisfying (4-4), 
there exist positive constants r] Q , p, and C ra d such that for 

1. any five points Z{ G IR? , 1 < i < 5, satisfying ( \j.3j ) with < rj < r] , and 

2. any function V : M 2 — > M 5 of the form l\4-o\) with Pert(0, 0) = and 
\\Pert(x,y)\\<pfor all \\(x,y)\\ < 1, 

the convex hull in M 5 of {V(0), V(z\), . . . , V(z$)} contains a ball of radius C ra &. 



Proof. Given the five points Zj G M 2 satisfying (|4.5| ), we define the 5x5 matrix 
N = N(ri,...,r B ,9 1 ,...,e B ): 



( T\ COs( 


Ox) 


r\ sin( 


Oi) 


1-2 
2 r l 


cos 2 ( 


9i) 


rj 


cos( 


<h) 


sin( 


9i) 


1-2 
2 r l 


sin 2 ( 


9i)\ 


r 2 cos( 


e 2 ) 


r 2 sin( 


9 2 ) 


1-2 
2 '2 


cos 2 ( 


92) 


r\ 


cos( 


62) 


sin( 


2 ) 


1-2 
2' 2 


sin 2 ( 


9 2 ) 


r 3 cos( 




r 3 sin( 


9 3 ) 


l r 2 
2 r 3 


cos 2 ( 


9s) 


r| 


cos( 


03) 


sin( 


9 3 ) 


1-2 
2 r 3 


sin 2 ( 


9 3 ) 


r 4 cos( 




r 4 sin( 


9a) 


l r 2 


cos 2 ( 




r\ 


cos( 


9a) 


sin( 


9,) 


1-2 
2'4 


sin 2 ( 


9,) 


\r 5 cos( 


e 5 ) 


r 5 sin( 


9b) 


1-2 
2'5 


cos 2 ( 


9 5 ) 


r\ 


cos( 


e s ) 


sin( 


9 5 ) 


1-2 
2'5 


sin 2 ( 


9s)J 



and a function ji: 

/i(ri, . . . , r 5 , 0i, . . . , B ) = — det jV(n, . . . , r 5 , 6>i, . . . , 5 ). 

Note that the constraints in ( |4.5| ) imply 

r 2 > C~ r \T 1 — t 2 > (1 — C + )r! 
r 5 > C~r 4 r 4 — r 5 > (1 — C + )r 4 
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so that 



M r i> r 2, r 3 , r 4 , r 5 , a x ,a u a 2 , a 3 , a 3 ) 
i 

r 1 r 2 rlr i r 5 (r 1 - r 2 )(r 4 - r 5 ) sin(a 2 - «i) sin 2 (a 3 - «i) sin(a 3 - a 2 ) 



4- (5! 

>^(C-) 10 (l-^sin (a2 - 
:= B(C _ ,C + ,a;i,a;2,a3) 



«i ) sin (a 3 — QfiJ sin(«3 — oti) 



Since determinants vary continuously with the matrix elements, there exists rjo 
such that if the Z; t satisfy (JO) with < r\ < r] Q , we have: 



M r i> r 2 ; r 3, r-4, r 5) 6»i, 6» 2 , 3 , #4, #5) > T^ 7 " 1 ' r2 ' r3 ' r4 ' r5 ' ai ' ai ' ° 2 ' " 3 ' 

> -jB{C^, C + , ai, a 2 , a 3 ) 

In particular, when < 77 < 770, we can bound /i from below independent of r iy 
8i, and rj. It follows that there exists p > such that if we perturb any matrix 
iV = N(ri, . . . , rs, 9\, . . . , 9$) by an error matrix E = (Ey) with \Eij\ < p, then 
det(iV + £) > ±det(N). 

Since V(0, 0) = 0, the simplex formed in M 5 from the six points V(0), V(z\), 
. . . , V{z§) has five-dimensional volume given by the determinant pp. 44]: 



1 , 

— det 

5 



(V( Zl )\ 



\V(zb)J 



> I_L d et (jV(n, . . . , r 5 , B x , . . . , 9 5 )) > ~B(C~ , C+, a 1} a 2 , a 3 ) 



whenever the Zi satisfy (|4.5|) with < rj < r] and \\Pert(x, y)\\ < p. Thus, the 
convex hull of these points has a guaranteed minimum volume. Since the positions 
of the points V(zi) are constrained in space, there is a constant C ra d > such 
that there will always be a small ball of radius C ra d contained somewhere within 
this convex hull. □ 



Lemma 4.10. Let C , C + , ai\, a 2 , a 3 , C ra d, Vo, ^nd p be as in Lemma Let 
Cpert > and < e < min (1, n p ) . Then, for 

\ ^ Pert / 

1. any five points Zi = (er, cos(#i), er, sin(#j)) in M 2 , 1 < i < 5, with n and 9i 
satisfying ( \4-3j) for some < r\ < r} , and 

2. any function V : M 2 —>■ M 5 of the form ( JJ> ) with 

\\Pert(x,y)\\ < C Pert \\(x,y)f for all \\(x,y)\\ < e, 
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the convex hull in M 5 o/{V(0), V(zi), . . . , V{z^)} contains a ball of radius C ra d€ 2 . 

Sketch of Proof. This proof is similar to that of Lemma |475| and follows from 
changing coordinate systems between the standard basis {ei, e 2 , e^, e±, e$} and 
the basis {eei, ee 2 , e 2 e3, e 2 e^, e 2 e^}. □ 



Proof of Theorem |4.8| . The hypotheses guarantee that the Eckmann-Ruelle 
linearization M(P) exists and is the unique best linearization as e k — > 0. We will 
show that there is a constant C (independent of k) such that \\M k — M(P) || < Ce k 
for large k. 

By Theorem [2.6| , there are e' > and C er > such that for every < e < e' 
we have W(M(P),P,e) < C er e 3 . Without loss of generality e\ < min{e', 1}. 
Then, 

W(M k , P, e k ) < W(M(P), P, e k ) < C er e\ for k = 1, 2, . . . 

For each k, let := — M(P) and let be the set of small tangent vectors 
at P: 

S k := {p, -PElR b :P 3 = vrfe) and \\p 3 - p\\ < ^} . 

For each Pj G Sk, \\Pj — P\\ < by ([2.1]) . As in the proof of Theorem [L2| , we 
see that 

II^U^II < 2C er e^ for each v £ Sk- 

Thus, we need only show that the convex hull Hull(Sk) contains a ball in M 5 . 
We will use Lemma [4. 10| to accomplish this. 



With respect to the canonical embedding basis at P, we have for Pj near P: 

Pi-P = (h a , h i2 , -h a , hah®, + Rem(ha, h i2 ,p) 

where Pi — p = (ha,h i2 ) G M 2 and Rem(ha, hi 2 ,p) is a vector-valued func- 
tion continuous in a neighborhood containing the closed ball of radius jh about 
(0,0). Define the function V(hi,h 2 ) as in Lemma [4.10|, with Pert(hi, h 2 ) = 



Rem(hi, h 2 ,p). Note that p — P = V(/iii, ^2) and we can take Cp erf = Crayior 
since ||Pem(/iji, ^2,p)|| < Cxayior \\(ha, hi 2 ) || 3 (see Chapter 2). It remains to 
construct the five points in M 2 that Lemma [4. 1 0| requires. 



By hypothesis, the attractor has three distinct fractal approach directions, 
and by Lemma [4.7| , these translate down to distinct fractal approach directions 
li, l 2 , 13 at p G iR 2 in the two-dimensional underlying attractor. These approach 
directions in M 2 are unit vectors, so we write li = (cos(aj), sin(aj)) for i — 1,2, 3. 
Since they are distinct approach directions, the angles a, are distinct and no two 
differ by a multiple of 7r. Moreover, we can arrange the li so that they satisfy (4.4). 
For each l i: there is a sequence of trajectory points in M 2 , q^ = p n (i,j),j = 1,2,..., 
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satisfying the definition of fractal approach direction for that Zj. Thus, there are 
constants < C~ < C + < 1 such that 



C~ < 



\Qij-P\ 



< C + for each i — 1, 2, 3 and j — 1, 2, 



(4.7) 



where the constants C and C + work for all three fractal approach directions. 
At this point, we have constants C~ , C + , a\, a 2 , and a^, and we obtain the 



constants C ra d, Vo an d P from Lemma 479. 



It will be convenient to write the points g^- in coordinates with origin at 
P- Qij ~ P — ( r ij cos(%), rij sin(%)). From condition (2) of the definition of 
fractal approach direction, we have y^Ipy — > U as j — ■> oo, which translates 
naturally to 0jj — > aj. We may now choose a number J large enough that 
for any i = 1,2,3 and j > J, — a.i\ < T) . Let K denote the least index 
k for which jr- < min (l, c £ - j and such that for each i = 1,2,3, we have 
maxdlgy : j > J} > ^. 

Now, let A; > if and consider the ball about p. Let gy, j > J, be the 
first point in the sequence converging to l\ for which \\qij — p\\ < Set 
z i — Qij ~P = ( r i cos (^i) ; r i sin(6 l 1 )). Since Hgi^—i) — p\\ > |f-, it follows from 
Q that n = ||<?ij > ^C - - Let z 2 = qiy+i) -p = (r 2 cos(6» 2 ), r 2 sin(6> 2 )) 
correspond to the next point in the sequence for l\. Thus, C~ < y < C + . Sim- 
ilarly, we choose points z± and z 5 for l 3 . We choose from l 2 in the same way 
we chose z\ and z$, i.e., to be the first point along the subsequence converging 
to l 2 within the radius of #-. We do not need a closer point along l 2 . Figure 4.2 
shows the configuration we have constructed. 

Let fi = and write Zi = (^fi cos(0j), j sin(#j)) for z = 1, ... ,5. By 
construction, the fi and Qi satisfy ( |4.5| ) with 77 = 770. Thus, we can apply Lemma 
4. 10| to conclude that the convex hull of the points V(0) = P, V(zi), . . . , V(z 5 ) 



contains a five- dimensional ball of radius %re|. Of course, each point V(zi 
corresponds to some vector Pj — P in Sfc, and it follows that Hull(Sk) contains 
this small ball. Finally, by Proposition [473], we conclude that 



AC, C 2 

\A k \\<^-^e k ioik>K. 

^rad 



□ 



We turn our attention now to formulating our convergence theorems in terms 
of least squares estimates rather than minimax estimates. As before, we need 
a way of comparing matrices, i.e., a way of measuring the error involved in the 
linearization process over a particular e-ball. In practice, this is done by summing 
the squares of the error terms. The techniques of least squares can then determine 
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the matrix which minimizes this sum of squared errors. Usually the data set is 
finite. However, we are interested in a convergence question which, by its very 
nature, requires an infinite data set. Thus, we must find a reasonable notion for 
least squares that allows us to consider infinite data. The arguments we give now 
are intended to motivate the definition that will follow. 

Imagine that we have exactly iV trajectory points Pq, . . . , Pjv-i m where 
F(Pi) = Pj+i, and we want to determine the best local linearization in the least 
squares sense over an e-neighborhood about P. We would then be looking for a 
matrix M which minimizes 

E m + i-F(p)-M(p l -p)\\ 2 

\\Pi-P\\<e 

where the sum is taken over only those trajectory points with ||Pj — P\\ < e. Any 
matrix which minimizes the quantity above will also minimize 

i N-l 

— £ \\P i+1 - F(P) - M(Pt -P)f XB(P,e)(Pi) 



i=0 



where N' is the number of data points e-close to P. This quantity represents 
the average of the squared errors. The characteristic function Xb{p,h) is used to 
selectively pick off just those values which are close to the base point P. For the 
infinite data case, we take the limit as N — ► oo, being careful to adjust the value 
of N' for each value of N: 



^ N-l 



where 



lim — £ \\P l+1 - F(P) - M(P t - P)f XB { P,e){Pi) (4. 

i=0 



N-l 

N' = E XB(P,e)( P i. 
i=0 



The number N' counts those trajectory points Pi (among the first iV iterates) 
which are e-close to P. We will rewrite (|4.8|) using the Birkhoff Ergodic The- 
orem [H|. In what follows, we let A be the attractor of the underlying sys- 
tem and /i its natural measure (an invariant, ergodic probability measure). Set 
I(P, e) = A H 7r _1 (i?(P, e)). For a matrix L, we define the differentiable function 

*l(x, P ) := Mf(x)) - 7r(/(p)) - L (ir(x) - n(p))f . 

For almost every po (with respect to /i), the limit (|4.8|) can be written: 



jfE^Lo^df {po),p)XB(P,e)(^(f (Po))) _ J A ^ L (x,p)xB(P,e)( 7r ( X ))M X ) 



lim 

jfTg^XBwW'iPo))) J A XB(P,e){< X ))M^ 



N 

1 



p(I(P,e))Ji(P,e) 



$ L (x,p)d/j(x) 
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Since the denominator fi(I(P,e)) is independent of the matrix L, any matrix 
which minimizes the limit in (|4.8j) also minimizes the quantity 



LS(L,P,e):=[ <S> L (x,n-\P))dfi(x)) . 

We say the matrix L\ is a better linearization of F over the e-ball about P 
than the matrix L 2 provided LS{L X , P, e) < LS(L 2 , P, e). 

For the Eckmann-Ruelle matrix M(P) of a reconstruction from M l into M m , 
it follows from the proof of the Local Linearization Theorem |2.2| that 

^M{p){x^~\P))dn{x) 



\\Remf(x — p) — M(P)Rem{x — p)|| 2 dfj,(x) 
<C 2 Taylor (l + \\M(P)\\) 2 f \x-p\ 2m+2 d^(x) 

J We) 

< C 2 Taylor (1 + \\M(P) ||) 2 ^ m+2 e 2m+2 /i (I(P, e)) 

since \x — p\ < C n \\Tr(x) — P\\ < C w e by (|2.1|) . Because the natural measure \i is 
a probability measure, 

LS{M{P),P,e) < Ce m+1 . (4.9) 

We can do better still if /i is a uniform measure on the unit interval [0, 1], that is 
/i is absolutely continuous with respect to Lebesgue measure and has a positive 
continuous density function !F{x). In this case, since I(P,e) is contained in the 
subinterval [p — C^e,p + C n e] by Q2.1Q , we have 

fi(I(P,e)) < [ P+C ^ dfx(x) < max iFix)] -2C^e. 

Jp-C^e 0<x<l 

Thus, we can improve fl4.9Q in this case to 

LS(M{P),P,e) < Ce m+ I. (4.10) 

The goal for the rest of this chapter is to prove that, in the case where \x is 
a uniform measure on [0,1], any sequence of matrices will converge to M(P) if 
each matrix in the sequence is a better linearization than M(P) over some small 
ball. Specifically, we prove: 

Theorem 4.11 (Least Squares Convergence, ]R l ^lR m ). Assume (Al) and 
(A2). Let A be a compact attractor for f with natural measure fi. Assume that \i 
is invariant, ergodic, and uniform, with positive density function T . Let P G M m 
be a limit point of the reconstructed attractor tt(A) for which (A3) and (A4) hold. 
Let {ek}^ =1 be a decreasing sequence of positive numbers, — > 0. Let {Mfc}^ 1 
be a sequence of matrices such that LS(Mk, P,€k) < LS(M(P),P,6k) for all k. 
Then M k -> M(P) as k -> oo. 
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Before giving the proof of |4.11| , we prove a few lemmas. In what follows, (■, ■) 
will denote the usual inner product of vectors in M m . We start with a definition. 

Definition 4.12. For a function f mapping the interval [-1,1] into M m , and for 
< e < 1, we define 

5(f,e) := min f f (h)) 2 dh. 

\\ri\\=l J-e 

Lemma 4.13. If C(h) = (h, \h 2 , . . . , ±h m ) for he [-1,1], then 5(C, 1) > 0. 



f 1 2 

Proof. Suppose 5(C, 1) = 0. Since the map 77 1— > / (r],C(h)) dh is contin- 
uous on the unit sphere in M m , there must be some unit vector v for which 
J (u,C(h)) 2 dh = 0. Hence, (u,C(h)) = for all h G [-1, 1]. However, v ^ 0, 

and so (z/, C{h)) = ^ — y/i n is a nontrivial polynomial, a contradiction. □ 



n=l ™ ! 



Lemma 4.14. Le£ C(h) be as above, and let Pert(h) be a continuous vector- 
valued function into M m such that \\Pert(h)\\ < Cp ert \h\ m+1 for h 6 [—1,1]. // 
V(h) := C(h) + Pert(h), then 5(V, e) > \5{C, l) e 2m+1 as e -> 0. 

Proof. Note that (17, V(/i)) 2 > (77, C(/i)) 2 + 2 (77, C(/V)) (77, Pert(h)}. Thus, 
5(V,e) > min f (77, C(/i)) 2 d/i + 2 min f (r],C{h)) (77, Pert(h)) dh 

\\ri\\=lj-e \\rj\\=lj- e 

> 5(C, e) - 2 min f (77, C(/i)) (77, Pert(h)} dh 

|MI=1 J-e 

We bound this minimum as follows, where e m is the unit vector (0, . . . , 0, 1): 

min f (r},C{h))(r],Pert(h)} dh < f (e m ,C(h)} (e m , Pert(h)) dh 

\\rj\\=lj-e ' J-e 

m! J-e 



(m + 1)! 



Next, we determine a lower bound for 5(C, e). Let E e denote the diagonal matrix 
E e = diag(e, e 2 , . . . , e m ), and note that 

H-Ee^H > e m IMI = e m for every unit vector rj. 
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Setting h = et and 7 = E t rj/ \\E t i]\\, we change variables in the integral defining 

1' (r],C(h)) 2 dh = e J ^(E e r),C(t)) 2 dt 

= e\\E iV \\ 2 f^,C{t)) 2 dt 
> e 2m+1 5(C, 1) 

and taking the minimum over all unit vectors 77, we obtain 

5(C,e) > 5{C,l)e 2m+1 . 

Since 5(C, 1) > 0, we have 

S(V, e) > 6{C, l)e 2m+1 - 2Cpert e 2m+2 > -5(C, l)e 2m+1 

(m + 1)\ 2 

for sufficiently small e. □ 



Lemma 4.15. There is a constant C\ > 0, independent of e, such that for all 
small e > 



p+e 



min / " (77, tt(x) -P) 2 dx> de 2m+1 . 

|M| = 1 Jp-e 



Proof. Recall that tc(x) — P = T(x — p) + Rem(x — p), where T(x — p) is the 
order-m Taylor polynomial: 

m 1 

T(x-p) = Y J -^ n \p)^-pT 



n=l 



nl 



and Rem(x — p) is the Taylor remainder term satisfying 

\\Rem{x-p)\\ < C Taylor \x - p\ m+1 . 

Let Q be the change-of-basis matrix such that with respect to the canonical 
embedding basis, Q~ l T{h) = (ji, \h 2 , . . . , ^r^ m ) p = C(h). Then, we can set 

7 := Q T r]/ Q T rj and V(x — p) := C(x — p) + Q~ 1 Rem{x — p), and we have 
(r), ir(x) - P) = (Q T r), C(x - p) + Q~ l Rem(x - p)) = \\Q T v\\ (7, V(x - p)) . 



In order to apply Lemma |4.14| to V, we note that because Q is a fixed nonsingu- 
lar matrix, \\Q~ 1 Rem(x — p)\\ < 1 1 <^ 1 1 1 CT ay ior\x — p\ m+1 ■ We will also need to 
examine Q T rj . Observe that Q T acts on the unit circle to produce an ellipse 
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with axes of length JWi, where cij > are the singular values of Q. Thus, for 
each unit vector 77, we have Q T i] > ^a min > where cr min = mincxj. Thus, 
setting h = x — p, we have for each unit vector 77 



[ P+€ (77, vr(x) - P> 2 dz = Q T V 2 f (7, V(h)) 2 dh > l-a mm 5(C, 1) 



_2m+l 



by Lemma |4.14j . The conclusion follows by taking the minimum over all unit 
vectors rj. □ 



Lemma 4.16. For each m x m matrix E, there is a unit vector r] e M m such 
that 

\\Ev\\ > 11^111(770^)1 

for all vectors v G M m . 



Proof. Let 770 be a unit eigenvector for E T E corresponding to the largest eigen- 
value a of E T E. It follows from the theory of the singular value decomposition 



that a 



EE 



E\\ [15, pp. 266]. Moreover, E E is self-adjoint, hence 



diagonalizable with an orthonormal basis of eigenvectors. For any vector v in 
M m , we can write v = (r] ,v) rj + w for some vector w perpendicular to i] . 
Since the various eigenspaces are orthogonal and invariant under E T E, we have 
(w,r] ) = (E T Ew,r] \ = (w,E T Er] \ = 0. Therefore, 



\Ev\ 



E T Ev,v 



( m ,v) 2 (E T E m , m ) + (E T Ew,w 
(r] ,v) 2 (0770,7/0) + \\Ew\\ 2 



> (vo,v) a 

The conclusion follows by taking square roots. 



□ 



Proof of Theorem The hypotheses guarantee that the Eckmann-Ruelle 

linearization M = M(P) exists and is the unique best linearization as — > 0. 
We prove that \\Mk — M\\ — > as k —> 00. For any k, we have 



I(P,e k ) 



\\{M k -M) (ir(x) - P)f dfi(x) 



< (LS(M k ,P } e k ) + LS(M } P,e k )) 2 
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where we have used the Holder Inequality for the last step. By hypothesis and 
(pl~T0l) , it follows that 



\\(M k - M) (tt(x) - P)\Ydn{x) < 4CV fc 



2 2m+3 



On the other hand, by Lemma L16 , there is a unit vector 7] such that 



\\{M k -M) (tt(x) - P)\\ 2 dfi(x) 



> \\M k -M\\ 2 [ ( Vo ,7r(x)-P) 2 dii(x). 



We need a lower bound in terms of e for the integral on the right. Note that 



P 



C rr 



p + 



Using the continuity 



the interval /(P, £&) contains the interval 

and positivity of the density function T and Lemma |L15| , we have for sufficiently 
small e k : 

(r) , ir(x) - P) 2 dfi(x) > \f{p) / ^ (rjo, ir(x) - P) 2 dx 



> \Hp)Ci 



2m+l 



for some constant C\ > 0, independent of e&. It follows that 
\\M k -M\\ 2 ^£lel^<4Chl^ 
and therefore \\M k — M\\ = 0(e k ) as e k — > 0, completing the proof. 



□ 
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Chapter 5 



Numerical Computations 



In this chapter, we describe numerical experiments illustrating the theory devel- 
oped in the previous chapters. As we shall see, we can compute numerically the 
Eckmann-Ruelle linearization in (noise-free) examples. The algorithm we use is 
similar to that outlined in 0, though the specific program was written to han- 
dle any function as a measurement function not just time-delay measurement 
functions. We briefly outline the algorithm here. 

Recall that the Eckmann-Ruelle procedure has three basic steps. In the first 
step, one reconstructs the attractor with a measurement function. In the second 
step, one determines a local linearization matrix at each point of the reconstructed 
trajectory. In the final step, one extracts the Lyapunov exponents from the 
linearization matrices obtained in the previous step. 

For our numerical experiments, we assume that we are given the underlying 
dynamical system / : M n —>■ M n and a measurement function 7r : M n —>■ M m 
for reconstructing the attractor. We iterate the function / a number of times 
to remove transients. Then, continuing to iterate /, we produce a trajectory 
{pi} € M n . Instead of storing the points pi, we apply the measurement function 
to each p^ and store the vectors P, = 7r(pj) G M m . We then sort our list of p by 
their first coordinates, keeping track of where each Pj ends up in the sorted list. 

The procedure for computing the local linearization matrix at a base point P 
is straightforward and based on least-squares methods. The goal is to determine 
the mxm matrix M which best satisfies M(P, — P)m P i+1 — F(P) for all Pj close 
to P. Given a neighborhood radius e, we quickly search our sorted list for those 
Pi whose first coordinate lies within a distance e of the first coordinate of P. Any 
other point Pj certainly lies outside the ball in M m of radius e about P. This first 
search often reduces the number of points we need to check carefully to less than 
10% of the trajectory. For each Pj in our shortened list, we construct the vectors 
Pj — P and P i+ x — P(P) and check their lengths. If both of these vectors have 
length less than or equal to e, then we keep them. Since the linearization we seek 
satisfies M (Pj — P) ~ Pj + i — P(P), we build matrices A and B for which the 
A;-th column of A is a vector Pj — P and the fc-th column of B is the corresponding 



43 



vector P i+ i — F(P). Then, our local linearization matrix M will be the m x m 
solution of the matrix equation MA = B. We solve this equation using a singular 
value decomposition routine taken from JTB . 



Once the local linearization matrices have been computed, we must com- 
pute the exponents. Recall from Chapter 3 that the Eckmann-Ruelle-Lyapunov 
(ERL) exponents of the reconstructed trajectory P ,Pi, . . . in M m are the values 
obtained by the limit 

h ER {Po^) ■= lim ~\n\\M n -. 1 M n - 2 ---MiM v\\ 

n— >oo fi 

for unit vectors v G M m , where = M(Pj) is the m x m local linearization 
matrix (i.e., the Eckmann-Ruelle linearization) at the point Pi in the trajectory. 
To evaluate this limit, we employ the treppen- iteration algorithm suggested in 
HH, |J. Given our sequence of linearization matrices, we use the QR matrix 
decomposition to find orthogonal matrices Qi and upper-triangular matrices R4 
(with non-negative diagonal elements) such that 

MiQ t ^ = QiRi for % = 0, 1, 2, . . . 

where we take Q-\ to be the m x m identity matrix. Then, we can write 

M„_iM n _ 2 • • • M x M = g n _i J R n -A-2 ■■■Ro- 

The orthogonal matrix Q n -i won't affect the matrix norm. Thus, 

h E ii(Po,v)= lim -In||i2 r ,_ 1 i2 r ,_2---i2ii2 i/||. 

The product R n -\ ■ ■ ■ Ro is upper-triangular, and its diagonal elements are the 
eigenvalues of the matrix, which we expect will grow like the Lyapunov exponents. 
As mentioned in Chapter 3, Theorem |3.2| justifies reading the exponents directly 
from the diagonal: 

J 71—1 

A fc = lirn - ^ In {(Ri)kk) for k = 1, . . . , m. 

n i=0 

Thus, we are able to compute the ERL exponents from our linearizations by using 
the diagonal elements from the QR decomposition. 
We shall now discuss several examples. 

In Chapter 2, we considered the doubling map on an interval reconstructed on 
the unit circle in M 2 where the underlying dynamics was f(p) = 2p (mod 2n) on 
[0, 2tt), and the measurement function was n(p) = (cos(p), sin(p)). We computed 
the Eckmann-Ruelle linearization in this case to be 

"<*.»>= L 2x(V + l)) (5 ' 1) 
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as a function of the points (x, y) on the unit circle. In numerical experiments, our 
routines to compute the best local linearization find good approximations to this 
matrix. For example, we can compute the linearization matrix at the fixed point 
(1,0) of the reconstructed dynamics F. Using a trajectory of 100,000 data points, 
the linearization matrix for a neighborhood of radius e = 0.05 is computed to be 



'3.999601 0.000000\ 
0.001048 1.999675 J 




M(1,0). 



Theorems [4.2| and |4 . 1 1| suggest that as the neighborhood radius decreases, the 



computed linearizations should converge in matrix norm to M(1,0). Figure |5] 
shows a graph of this convergence as the radius shrinks to zero. To generate 
this graph, we first compute and store a reconstructed trajectory of 100,000 data 
points on the unit circle. Then, for each radius e, we find those points in the 
trajectory within the specified radius of our base point and use them to compute 
the linearization matrix. 

Likewise, the local linearizations converge at other points on the unit circle. 
For a neighborhood radius of e = 0.05, the linearization matrix computed at 
(^,^) « (.707, .707) is 



'1.414398 -1.413499\ _ / ^2 -^2> 
^2.828112 2.828112 ) ~ \2^2 2y/2 j 



M 




Figure |5]2| shows the convergence with radius of the linearization matrices at this 
point. 

For a given neighborhood radius, we can determine at each point on the circle 
the error between the computed linearization and the Eckmann-Ruelle lineariza- 



tion. Figure |5^ shows this graph for several different radii. 

With close agreement at each point in the interval between the numerically- 
determined linearization and the Eckmann-Ruelle linearization, we expect that 
the computed Lyapunov exponents will match those given by the Lyapunov ex- 
ponent formula of Theorem |3~3| . To test this, we computed the Lyapunov expo- 
nents for this example using 100,000 data points on the circle. We found values of 
0.693173 « m(2) and 1.386244 « 2 ln(2), as predicted. The first value is the true 
Lyapunov exponent of the doubling map, while the second value is a spurious 
exponent. Notice that the largest computed value is not a Lyapunov exponent 
of the underlying dynamical system. 

The Lyapunov exponent formula of Theorem |3.3| holds for any generic map and 
specifically for the one-to-one maps with linearly independent derivative vectors. 
Since the formula itself is independent of the measurement function n, we expect 
to compute the same exponents when different measurement functions are used. 
Changing the measurement function will likely change the speed of convergence 
somewhat, but we still expect the results to be close. With this in mind, we 
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Figure 5.1: Graph of the difference (in matrix norm) between the computed 
linearization and the Eckmann-Ruelle linearization at the fixed point (1,0) for 
the doubling map f(p) = 2p(mod 2tt) reconstructed on the unit circle. The 
calculation used 100,000 data points. Note the downward trend of the graph, 
indicating the convergence of the linearizations to the Eckmann-Ruelle matrix. 
The spikes in the graph appear to be numerical artifacts. For radii below about 
0.01, very few data points lie sufficiently close to the base point to use in the 
calculation. 
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Figure 5.2: Graph of the difference (in matrix norm) between the computed 
linearization and the Eckmann-Ruelle linearization at the point ^) for the 
doubling map reconstructed on the unit circle. The calculation used 100,000 data 
points. Note the downward trend of the graph, indicating the convergence of the 
linearizations to the Eckmann-Ruelle matrix. The spikes in the graph appear 
to be numerical artifacts. For radii below about 0.01, very few data points lie 
sufficiently close to the base point to use in the calculation. 
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Figure 5.3: These graphs show, for different radii, the error in matrix norm 
between the computed linearization and the Eckmann-Ruelle linearization for 
the doubling map at each point on the unit circle (expressed as an angle). As the 
radius decreases, the error decreases uniformly. The graphs were created using a 
single trajectory of 200,000 data points on the unit circle. 
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Figure 5.4: The curve ir(p) = (cos(27rp), sin(27rp), \e v * p ) for p e [0, 1]. 



performed a Lyapunov exponent calculation for the doubling map reconstructed 
on an ellipse using the measurement function 7i(p) = (2 cos(p), sin(p) — |cos(p)) 
and found exponents of 0.693151 « ln(2) and 1.386267 « 2 ln(2). 

Next, we examine the case of a reconstruction from iR 1 into M 3 . We use 
the logistic map f(p) = 4p(l — p) on [0, 1] reconstructed onto a curve by the 
measurement function ir(p) = (cos(27rp), sin(27rp), |e 1_ 2P). Figure |5.4| shows this 
curve in M 3 . 

In this case, the formula for the Eckmann-Ruelle linearization is more com- 
plicated than (|5.1|) , but for any individual base point, it is easily computed 
from the definition given in Chapter 2. For example, at the point 7r (0.125) ~ 
(0.707,0.707, 1.277) on the curve in M 3 , the Eckmann-Ruelle linearization is 



M(tt(0.125)) 



/-0.395832 2.313014 -18.147168 N 

-6.279704 -3.080283 65.170183 
V-0.000276 -0.098277 1.550445 



and, computing with neighborhood radius e = 0.05, the local linearization matrix 
is 

/-0.391945 2.313262 -18.187396^ 

-6.273483 -3.079227 65.103345 
V-0.000214 -0.098268 1.549787 

Figure |575| shows the error in the linearizations over the whole curve. To compute 
the Lyapunov exponents for this example, we used 100,000 data points on the 
curve. We found values 0.6884401 « ln(2), 1.392024 « 21n(2), and 2.056767 w 
3 ln(2). This is consistent with Theorem |3.3| since the Lyapunov exponent for the 
logistic map is ln(2). Note that the largest computed exponent is spurious. 

We also computed the Lyapunov exponents of several other reconstructions 
of the logistic map in M 3 . The results are listed in Table 5.1 . 

In addition to the somewhat arbitrary measurement functions we have used 
thus far, we can also construct time-delay embeddings in M 3 of the logistic map. 
At each iteration, we record some quantity based on the current state of the 
system. In an actual experiment, this choice would be realized by the appara- 
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Radius = 0.05 




Figure 5.5: These graphs show, for different radii, the error in matrix norm 
between the computed linearization and the Eckmann-Ruelle linearization of the 
logistic map at each point on the curve in M 3 . As the radius decreases, the error 
decreases uniformly. The graphs were created using a single trajectory of 200,000 
data points on the curve in M 3 . 



Table 5.1: The Eckmann-Ruelle-Lyapunov exponents of several reconstructions 
of the logistic map f(p) = 4p(l —p) into M 3 . Each computation involved 100,000 
data points. The computed exponents are roughly ln(2), 21n(2), and 31n(2). 



Measurement Function n{p) 


Computed ERL Exponents 


(p,p 2 ,p 3 ) 
(p, ^/p - 2p 3 , |pcos(2.57rp)) 
(cos(27rp),p + sm(2irp) , p + p sm(5np 2 )) 


0.664767 1.411537 2.082149 
0.678497 1.391379 2.075309 
0.687375 1.385631 2.066531 
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Table 5.2: The Eckmann-Ruelle-Lyapunov exponents of several time-delay recon- 
structions of the logistic map f(p) = 4p(l — p) into IB?. For each iteration, we 
recorded a scalar determined by the current state p of the system. Each com- 
putation used 100,000 data points. The computed exponents are roughly ln(2), 
21n(2), and 3m(2). 



Recorded Quantity 


Computed ERL Exponents 


p 

cos(2p) — p 
ln(l +p) 


0.691530 1.386031 2.083390 
0.692866 1.386281 2.080733 
0.691398 1.386725 2.079474 



tus used to measure and record the data. The results of computing Lyapunov 
exponents for several time-delay examples are given in Table |5.2| . 

The Eckmann-Ruelle linearization described in ( |5.1|) may be misleadingly sim- 
ple. It is possible to compute an explicit formula for the Eckmann-Ruelle lin- 
earization for the case of the logistic map f(p) = Ap(l —p) on [0, 1] reconstructed 
onto the unit circle by the measurement function ir(p) = (cos(27rp), sin(27rp)). 
One then finds that some of the component functions of that matrix can have 
derivatives as large as 200 in absolute value. See Figure |5~6| . 

We move on to discuss the case of two-dimensional underlying dynamics. 
Here, we base our experiments on the Henon map f(x, y) = (1.4 — x 2 + 0.3y, x). 
We begin by considering the measurement function n(x,y) = (x,y,x 2 ,y 2 ,xy). 
For this measurement function, we can compute the Eckmann-Ruelle linearization 
explicitly: 



/ 



M(ir(x,y)) 




1 

1.2x2/ — 8x 3 


1.4 + 3x 2 



0.3 -1 


0.84 + 0.6x 2 -2.8 + 6x 2 - 0.6y 
1 
-3a; 






.09 








-1.2a; 


0.3 



where for convenience, we refer to M using coordinates in the underlying space 
M 2 instead of coordinates in the reconstruction space M 5 . We observe the same 
phenomena when the underlying dynamics are two-dimensional as we did in the 
previous cases. For example, at the point 7r(l. 555478, 0.398567) G M 5 on the 
reconstructed Henon attractor, the linearization matrix for a neighborhood of 
radius e = 0.01 is computed to be 



0.000000 
1.000000 
-29.365725 
0.000000 
8.657187 



0.300000 
0.000000 
2.290703 
0.000000 
-0.001385 



-1.000000 
0.000000 

11.478416 
1.000000 

-4.666094 



0.000000 
0.000000 
0.090305 
0.000000 
0.000326 



0.000000 \ 
0.000000 
-1.866085 
0.000000 
0.300723 J 
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Figure 5.6: A graph of the first coordinate of the Eckmann-Ruelle matrix for the 
logistic map /(p) = 4p(l — p) on [0,1] reconstructed on the unit circle in M 2 by 
the measurement function ir(p) = (cos(2iTp), sin(27rp)). Note the large derivatives 
near the ends of the interval. 



At another point, 7r(— 1.741541, 1.753985), on the reconstructed Henon attractor 
in M 5 , the local linearization for e = 0.01 is computed to be 

/ 0.000000 0.300000 -1.000000 0.000000 0.000000\ 

1.000000 0.000000 0.000000 0.000000 0.000000 

38.590278 2.658478 14.345224 0.090317 2.089740 

0.000000 0.000000 1.000000 0.000000 0.000000 

V10.498785 -0.000234 5.224592 0.000067 0.300000/ 

Both of these examples are in good agreement with the general formula 
above. As before, we can graph the convergence of the local linearizations to the 
Eckmann-Ruelle linearization as the neighborhood radius shrinks to zero. See 
Figure pVT[ Using 300, 000 data points, we computed the Lyapunov exponents 
for this example. Recall that the true Lyapunov exponents of the Henon map 
are approximately A = 0.42 and p, = —1.62. We found values of 0.886200 w 2A, 
0.422459 « A, -1.195450 « A + p, -1.636780 « p, and -3.183227 « 2p. These 



values are consistent with the Lyapunov exponent formula given in Theorem |3T6 . 

We also computed the Lyapunov exponents from time-delay reconstructions 
of the Henon system. Table |5.3| shows the results. As one can see in Table |5.3|, 



we do not find all of the exponents predicted by the Lyapunov exponent formula 
in Theorem |3.tj . Three of the exponents match nicely with the predicted formula, 
namely A ~ 0.42, 2A ~ 0.84, and p ~ —1.62, but the other two vary somewhat 
from their expected values. With more data and smaller neighborhood radii in 
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\ 

0.12 




0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 



Radius 



Figure 5.7: Graph of the difference (in matrix norm) between the computed 
linearization and the Eckmann-Ruelle linearization at the point (-1.741541, 
1.753985) for the Henon map reconstructed in M 5 . The calculation used 1,000,000 
data points. The graph does not show radii larger than about 0.04 because this 
base point lies about 0.045 units away from a significant bend in the Henon 
attractor. 
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Table 5.3: The Eckmann-Ruelle-Lyapunov exponents of several time-delay recon- 
structions of the Henon map into M 5 . For each iteration, we recorded a scalar 
determined by the current state (x, y) of the system. Each computation involved 
300,000 data points. Three of the computed exponents match the Lyapunov ex- 
ponent formula: 0.84 « 2A, 0.42 m A, and —1.62 « fi. The other exponents will 
probably converge for longer time series. 



Recorded Quantity 


Computed ERL Exponents 


X 


0. 


338426 


0.41* 


3513 


-1.070009 


-1.619577 


-2.288856 


y-xy 


0. 


339063 


0.41* 


3518 


-1.064124 


-1.619594 


-2.426907 


\ (x 2 + y 2 ) 


0. 


337956 


0.41* 


3331 


-0.846145 


-1.621547 


-2.404775 


arctan 


0. 


340121 


0.41* 


3186 


-1.062045 


-1.618537 


-2.226276 


cos(x) + sin(y) 


0. 


339557 


0.41* 


3119 


-1.019653 


-1.621305 


-2.378717 


cos(y - 1) - fx 2 


0. 


339295 


0.41* 


3780 


-1.095620 


-1.621523 


-2.624901 



the computations, those last two exponents would likely converge to their proper 
values. Note that the true Lyapunov exponents appear in each of these examples, 
suggesting that they converge to their correct values rather quickly. This is 
not surprising if one looks back at the proofs of Theorem |2.6| in Chapter 2 and 
Theorem |3.6] in Chapter 3. Any matrix that does not map the first-order terms 
of the Taylor series correctly will have 0(||AP||) error, instead of the 0(||AP|| 3 ) 
error of the Eckmann-Ruelle linearization. Fortunately, algorithms to compute 
local linearizations will instead find matrices with error 0(||AP|| 2 ) or better, and 
these matrices will map the first-order Taylor terms correctly. This ensures that 
the true Lyapunov exponents will be among the first of the computed exponents 
to converge. 



54 



Appendix A 
Polynomial Results 



In this appendix, we derive the facts about polynomials that we need in the text. 

Proposition A.l. Let p(x) be a polynomial of degree d > 0, and suppose there 
are infinitely many values Xi — > for which \p(xi)\ < C\xi\ q for some q > d. 
Then p(x) = for all x. 



Proof. Write p(x) = adX d + h a^x + a . By the continuity of p(x), 

< |oo| = |p(0)| = lim \p(xi)\ < lim C\xi\ q = 

since q > d > 0. Thus, a = 0. 

For induction, suppose that oq = a± = ■ ■ ■ = a^-i = where k < d < q. 
Write p(x) = x k {adX d ~ k + • • • + a^+ix + a k ) := x k qk(x). Then, by the continuity 
of q k (x), 

< \a k \ = \q k (0)\ = lim \q k ( Xi )\ = lim < lim ^! = ^ m C \x i \^ k = 

since q — k > q — d > 0. Thus, a k = 0, and the lemma follows from induction on 
k. □ 

We would like to prove a similar result for polynomials with more variables. 
Unfortunately, the result is not true, even for two variables, without extra hy- 
potheses. For example, if p(x,y) = x — y, then p(xi,Xi) = for any sequence 
Xi — > 0. Thus, p(x,y) need not be zero, even though there is a sequence going 
to (0,0) satisfying \p(xi,yi) \ < C\\(xi,yi)\\ g for every q > 0. The added hypothe- 
ses in the next lemma guarantee that the polynomial must be small in many 
directions, not just one unlucky choice. 

Recall the definition of an approach direction from our discussion in Chapter 
2 of the case where the underlying dynamics is two-dimensional. 
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Definition A. 2. Let I be a unit vector in IR m . A subset S of IR m has the ap- 
proach direction I at the base point P G M m provided there is a sequence 
{QfcjfcLi from S such that: 

1- Qk — > P (is k — > oo ; and 

2. ii^q^i I as k — ► oo, where £Q k := Q k — P. 

We call a collection of approach directions at P distinct provided that no two 
are the same and no two are reflections through the origin. Multiple approach 
directions are not required to be linearly independent. 

Proposition A. 3. Let p(x,y) be a polynomial of (total) degree d (that is, each 
monomial in p has total degree at most d). Assume that there is a sequence of 
points (x k , Vk) (0, 0) in M 2 such that \p(xk, y k ) | < C\\ (x k , y k ) \\ q for some q > d. 
If the set {(x k ,y k ) : k > 0} has d + 1 distinct approach directions at (0,0), then 
p(x, y) = for all (x, y) G M 2 . 

Proof. Write p(x,y) = T,j=oPj{x,y), where Pj{x,y) = Y,i=o a jeX e y J ~ e is a 
polynomial of degree j with only monomials of total degree j. We will prove 
inductively that each polynomial pi(x, y) is identically zero for I — 0,1, ... ,d. 
We begin with the I = case. Note that 

< |a 00 | = b(0,0)| = lim \p{x k ,y k )\ < lim C\\ (x k , y k ) f = 

since q > d > 0. Thus, po(x, y) = aoo = 0, and so p has no constant term. 

Without loss of generality, we may now assume that d > 1. 

Next, we show that p±(x, y) = a n x + a w y = by proving that an = a w = 0. 
This is the / = 1 case and it will demonstrate the basic idea for the induction that 
follows. For an approach direction (a, (3), let {(x k ,y k )} be a subsequence such 
that (x k ,y k )/\\(x k ,y k )\\ -> {a, (3). For convenience, set n k := \\{x k ,y k )\\. Then, 

< \a±ia + a w /3\ = lim 

fc^oo 

, v ( \p(x k ,y k )\ , i \x k \ e \y k \ j - e \ 

< hmsup h > > a,-J 

fc-oo y n k n k J 

< limsup [ + ^ kiiPfcP" 1 — + \aje\\x k \ e \y k \^ e " 1 — J 

fe^oo y n k - =2 n k J=2 e= o n k J 

= 



Xk . y k 

an h aio — 

n k n k 



lim 

fc^oo 



\pi{xk,yk) 

n k 
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because q > d > 1, \xk\/nk — > «, li^l/n^ — > /3, and jy^] — > 0. Therefore, 
a n a + a w j3 = for each approach vector (a, (3). When a ^ 0, we can divide by 
a to get: 



\ a . 







We must now consider several cases. 

CASE 1: Assume that either d > 2, or d = 1 but there are 2 distinct approach 
vectors with a ^ 0. In this case, there are at least two approach vectors (ct,f3) 
with a / 0. If we write f(t) = aiot + an, then for each such vector, we have 
/ = 0. As there are at least two such values and f(t) is linear, we must 
have f(t) identically zero. Hence, an = a w = and therefore pi(x,y) = for all 
(x,y). 

CASE 2: Assume that d = 1 and that one of the 2 approach directions is 
(0, ±1). In this case, we see directly that |aio| = |an0 + aio(±l)| = 0. Then, for 
the approach vector (a, (3) with a ^ 0, we have = janct + • f3\ = \a,ua\. We 
conclude that a n = as well. Therefore pi(x,y) = for all (x,y). 

In all cases, pi(x,y) = as desired. The main induction step that follows is 
similar. 

Suppose we have shown that po(x,y) = pi(x,y) = ■■■ = pi-i(x,y) = are 
all identically zero for some / < d. We want to show that pi(x,y) = is also 
identically zero. For any approach direction (a, (3) (using the same notation as 
above): 



< 



E a le® e (3 
e=0 



e ol—e 



lim 

fc^oo 



^0 \ n k) \ n k) 



= lim 

k— >oo 



\Pi{xk,Vk) 



fc^oo 



<\imsup[Cnl- l + E E T : 



3=1+1 
d 3 



ni 



k— >oo 



j=Z+l e=0 



where the terms Tt look like 



e-piL'; \j—e—p2 



I - I \ PI / I * I \ P2 



for some appropriate choice of pi + £>2 = £ < j- The lim sup above is zero because 

1 

q-l>q-d>0, —(x k ,y k )-*(a,P), and \x k \,\y k -> 
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together imply that for each j and e, the terms T^ e — > as k —>■ oo. Thus, 



53 aiea e (3' 







for each (a, (3). 



e=0 



For those approach vectors with a ^ 0, we can divide through by a 1 : 




0. 



(A.l) 



We now must consider a few cases. 

CASE 1: Assume that either d > I, or d = I but there are d + 1 distinct 
approach vectors with a ^0. In this case, there are at least 1 + 1 approach vectors 
(a, j3) with cr ^ 0. If we write f(t) = Yj l e =o a iet l ~ e then for each such approach 



vector, we have / (f J = by flOp . As there are at least I + 1 such values and 

deg (/(£)) < I, we must have f(t) identically zero. Hence, aio = ■ ■ ■ = an = and 
therefore pi(x,y) = for all (x,y). 

CASE 2: Assume that d = I, and that one of the I + 1 approach directions is 
(0, ±1). In this case, we see directly that 



Thus, a/o = and the polynomial f(t) defined as in the previous case is actually 
of degree at most I — 1. Then, for the d = I approach vectors (a, (3) with 
we have / ( £j = 0. We conclude, as before, that f(t) is identically zero. Hence, 
a l0 = ■■ ■ = a a =0 and pi(x, y) = for all (x, y). 

In all cases, pi(x,y) = as desired. This completes the induction step and 
the proof. □ 




a l0 \ = 53 a *eO e (-l) 



l-e 







e=0 
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Appendix B 



Lyapunov Exponents of Upper Triangular 
Matrix Products 



In this appendix, we prove Theorem 372 that when a sequence of upper-triangular 



m x m matrices are multiplied together, we can find vectors whose lengths grow 
at the same rate as the diagonal elements of the matrix product. This result 
can be found elsewhere in the literature, for example in [I7J , though it usually is 



presented in a more general context. This presentation of the theorem roughly 
follows [I7j but can also be considered a distillation of a piece of Oseledec's 



original proof. 

Lemma B.l. For k = 1,2,..., let A k be an m x m matrix written in block 
structure A k = ( a fe |* ), where a k e JR, b k e #? lx ( m - 1 ) j and B k e ]R(rn~i)x(m-i) _ 
Set Sq = I m , the m x m identity matrix, and for n = 1,2,... define S n = 
A n - ■ ■ A\ = ( S q j^j) using the same block structure. Then 

, _ ^ b k T k _i 

t n ci n • • • a,\ y . 

k=1 a-k ■ ■ ■ o-i 



i b 

Proof. The proof is by induction on n. For n — 1, t\ — a\ y^ — = b\. In 



general, 



k=l 



«1 



t n +l — a n+\tn + b n +iT n 

E n b k T k 
r a n +i ■■■ai 

k=l a k' ' ' a l dn+l • • • Ol 



' L b k T k _i 



A.- 



1 a k - ■ ■ a\ 



□ 
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Definition 3.1. 

7 provided 



A sequence of real numbers {?"•&} has (geometric) growth rate 
In IrjJ 



lim 

k— >oo 



7- 



Theorem 3.2. For = 1, 2 ... 7 let Af. be m x m upper triangular matrices, and 
define S n = A n ■ ■ ■ Ai. Assume the magnitudes of the entries of A^ are bounded 
independent of k, and that the diagonal entries of S n have growth rates 7i, . . . ,7 m 
as n — > oo. TTien i/jere exists vectors v±, . . . ,v m such that for each i = 1, . . . , m, 
\\S n Vi\\ has growth rate 7$. 

Proof. The proof is by induction on m. Notice that when multiplying two 
upper triangular matrices, the lower i rows of the product are independent of the 
entries of the two matrices above the lower i rows. So, the induction on m will 
start at the lower right and work upward. 

Case m — 1 is clear; the vector v\ = (1) works. For the general case, assume 
that the theorem holds for all such sequences of (m — 1) x (m — 1) matrices. We 
use the partitioned notation for A^ and S n from Lemma p.l| . By the induction 
hypothesis, we will assume that by denoting by v 2 , . . . ,v m the columns of the 
matrix 

A v 2 3 ... v 2m \ 



V, 



m—1 



V3r, 



V 1 / 

the sequence ||T n z)j|| has growth rate 7^ for % = 2, . . . ,m. (Recall from Lemma 



B.l| that T n is the lower right (m — 1) x (m — 1) submatrix of S n .) In particular, 
for each e > there exists a constant K such that 



for i = 1, . . . , m, and furthermore, by assumption 

1 



K 



e (7i-e)n < Sn = a n ---a 1 < Ke ( ^ 1+e)r 



(b.i; 



(B.2) 



We will add a top row to V m -\ to get 



V„. 



(I V12 ... Vi m \ 
1 • • • V 2m 



V 



1 / 



such that using the new columns v\, . . . , v m of V m , the sequence HS^H has growth 
rate 7« for i — 1, . . . , m. In fact, this is already clear for z = 1; we just have to 
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make sure none of the other growth rates were changed by adding the top row of 
V m . 

Now we give the definitions of the new entries in the top row, v i 2 , . . . , v ± m . If 
71 < 74, set vu = 0. If 71 > 7*, define 



Vli 



E 



b k T k _iVi 



k=l a k ' ' ' a l 



(B-3) 



This series converges by comparison to a geometric series, because for each e > 0, 
the magnitude of the kth term is bounded above by a constant (independent of 
k) times e (^- e ) k = e( 7 *~ 7l+2e ) fc , using equations (|B.1|) and (|B.2|) . Here we used 
the assumption that the magnitude of the entries of b k from A k obey a bound 
independent of k. 

To finish, we need to show that for each i = 2, . . . , m, the growth rate of 
llS^uJ is 7j. There are two cases. First, we consider i such that 71 < 7^. 



S n Vi 



tnVi 

. T n Vi . 



_ s^n b k T k _ 1 v i ' 



T n Vi 



(B.4) 



where we have used Lemma B.l to rewrite t n . Equations ( |B.1|) and ( B.2 ) imply 
that for each e > there exist constants independent of n such that 



EbkTk-iVi 


k=1 &k ■ ■ ■ Oi 



< constant ■ e (71+e) " V - 7 - 

An 

< constant ■ e (7l+e)n 



fc=i 

e ( 7l - 7l +2e)(n+l) _ 1 



< constant ■ e (7l+e)n 

< constant ■ e (7i+3e)n 



e 7i-7i+2e _ X 
e (7 l -7i+2e)n g(7i-7l+2e) 



e 7i-7l+2e _ X 



Since both entries of S n Vi have 7^ as an upper bound for the lim sup of the growth 
factor as n — > 00 (using the above inequality and equation ( |B.1[ )), and since the 
lower entry has 7^ as a lower bound for the lim inf of the growth factor (again 
from (|B.1| )), the growth factor limit exists and is 7^. 
Second, we treat the case where 71 > 7^. 



S n V\i + t n Vi \ 
K T n Vi J 

-a n ■ ■ ■ a x 2^ = i ak ... ai +a n - ■■a 1 l^ k=1 



a k ---ai 



a l J^k^n+l 
T n Vi 



T n Vi 

b k T k _ 1 v i - 
a k ---ai 
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We can bound the top entry as follows. For each e > there exist constants 
independent of n such that 



k=n+l ak 



hTk-iVi 



ai 



oo ( 7i +e)fc 

< constant ■ e (71+e)n V -. -r 

fc=n+l e 



= constant ■ e 



(7i+e)n 



= (7 l -7l+2e)(n+l) 



1 — e 7 l ~7l+2e 



< constant • e (7l+3e)n . 



Just as in the first case, both entries of S n Vi have 7$ as an upper bound for the 
growth factor, and the lower entry has 7$ as a lower bound for the growth factor. 
Thus, the growth factor is 7*. □ 
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